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research,  and  these  reports  are  consolidated  into  this  annual  report. 
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PREFACE 


This  volume  is  part  of  a  16-volume  set  that  summarizes  the  research  accomplishments  of 
faculty,  graduate  student,  and  high  school  participants  in  the  1992  AFOSR  Summer  Research 
Program.  The  current  volume,  Volume  11  of  16,  presents  the  final  research  reports  of  graduate 
student  (GSRP)  participants  at  Arnold  Engineering  Development  Center,  Civil  Engineering 
Laboratory,  Frank  J.  Seiler  Research  Laboratory,  and  Wilford  Hall  Medical  Center. 

Reports  presented  herein  are  arranged  alphabetically  by  author  and  are  numbered 
consecutively  -  e.g.,  1-1,  1-2,  1-3;  2-1,  2-2,  2-3.  Research  reports  in  the  16-volume  set  are 
organized  as  follows: 

VOLUME  TITLE 

1  Program  Management  Report 

2  Summer  Faculty  Research  Program  Reports:  Armstrong  Laboratory 

3  Summer  Faculty  Research  Program  Reports;  Phillips  Laboratory 

4  Summer  Faculty  Research  Program  Reports:  Rome  Laboratory 

5A  Summer  Faculty  Research  Program  Reports;  Wright  Laboratory  (part  one) 

SB  Summer  Faculty  Research  Program  Reports:  Wright  Laboratory  (part  two) 

6  Summer  Faculty  Research  Program  Reports:  Arnold  Engineering  Development  Center;  Civil 
Engineering  Laboratory;  Frank  J.  Seiler  Research  Laboratory;  Wilford  Hall  Medical  Center 

7  Graduate  Student  Research  Program  Reports:  Armstrong  Laboratory 

8  Graduate  Student  Research  Program  Reports;  Phillips  Laboratory 

9  Graduate  Student  Research  Program  Reports:  Rome  Laboratory 

10  Graduate  Student  Research  Program  Reports;  Wright  Laboratory 

11  Graduate  Student  Research  Program  Reports:  Arnold  Engineering  Development  Center;  Civil 
Engineering  Laboratory;  Frank  J.  Seiler  Research  Laboratory;  Wilford  Hall  Medical  Center 

12  High  School  Apprenticeship  Program  Reports:  Armstrong  Laboratory 

13  High  School  Apprenticeship  Program  Reports:  Phillips  Laboratory 

14  High  School  Apprenticeship  Program  Reports:  Rome  Laboratory 

15  High  School  Apprenticeship  Program  Reports:  Wright  Laboratory 

16  High  School  Apprenticeship  Program  Reports:  Arnold  Engineering  Development  Center;  Civil 
Engineering  Laboratory 
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ibstract 


Experience  using  a  model-based  approach  to  develop  an  83  processor  parallel  instrumen¬ 
tation  system  for  turbine  engine  aeromechanic  stress  analysis  is  described.  The  approach 
includes  using  a  graphics  based  editor  to  describe  the  structure  of  the  desired  signal  flow 
graph  as  well  as  the  target  hardware  architecture.  Program  synthesis  techniques  are  used 
to  automatically  transform  these  models  into  an  executable  system. 
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AN  OVERVIEW  OF  THE  BEHAVIOR  OF  ALUMINUM 
IN  SOLID  PROPELLANT  ROCKET  MOTORS 


Robert  J.  Geierman 
Graduate  Research  Assistant 
Department  of  Aeronautical  Engineering 
University  of  Tennessee  Space  Institute 

ADatract 

An  extensive  literature  search  was  conducted  in  order  to  provide 
an  adequate  understanding  of  the  processes  which  aluminum  undergoes  at 
several  locations  in  a  solid  rocket  motor.  This  paper  describes  the 
phenomena  which  occur  at  the  propellant  surface,  Che  combustion  chamber, 
the  rocket  motor  nozzle,  and  the  exhaust  plume.  These  descriptions 
include  a  discussion  of  previous  models  and  experiments  that  have  been 
conducted.  Although  several  of  these  models  make  very  accurate 
performance  predictions,  much  of  their  basis  rests  on  emperical  data 
instead  of  analytical  models.  Due  to  this  fact,  some  of  the  previous 
models  may  have  several  shortcomings (in  the  analytical  sense).  Some  of 
these  shortcomings  include:  -1)  the  lack  of  an  adequate  analytical 
agglomeration  model,  -2)  the  neglection  of  agglomerate  radiation  heat 
transfer  to  Che  propellant  surface,  -3>  the  neglection  of  particle 
collisions  and  fragmentations,  -4)  no  predictions  of  slag  accumulation 
versus  nozzle  geometry,  and  -5)  the  lack  of  an  accurate  description  of 
the  particle  distribution  and  nozzle  ablation  at  the  nozzle  throat.  An 
effort  to  remedy  these  shortcomings  will  be  presented  in  a  thesis  at  a 
later  date. 
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Introduction 


Metal  additives  have  long  been  used  in  solid  rocket  propellants  in 
order  to  improve  acoustic  stability,  control  burn  rate  and,  most 
importantly,  to  improve  specific  impulse.  Due  to  its  relatively  low 
cost,  availability,  and  high  combustion  temperature,  aluminum  is,  by 
far,  the  most  often  used  metal  additive. 

Although  aluminum  has  been  used  in  rocket  motors  for  many  years, 
its  combustion  processes  are  still  not  completely  understood.  An 
adequate  understanding  of  the  combustion  of  aluminum  is  necessary  to 
accurately:  -1)  predict  plume  signatures  and  radiation  heat 

transfer (due  to  the  aluminum  oxide  particle  radiation  in  the  plume), 

-2)  predict  combustion  stability  due  to  acoustic  damping  of  aluminum 
oxide  particles,  and  -3)  compute  the  specific  impulse  of  the  rocket 
motor  by  properly  incorporating  the  two-phase  flow  losses  of  the 
aluminum  oxide  particles. 

This  paper  attempts  to  describe  the  physical  processes  which  occur 
throughout  the  aluminized  solid  propellant  rocket  motor.  This  is 
accomplished  through  the  use  of  an  extensive  literature  search  which 
describes  the  physical  phenomena  for  four  locations  throughout  the 
motor.  These  locations  include:  -1)  the  propellant  surface,  -2)  the 
area  inside  the  combustion  chamber,  -3)  the  rocket  motor  nozzle,  and 
-4)  the  rocket  motor  exhaust  pl'ime.  Later,  these  processes  will  be 
described  by  partial  diff- rential  equations,  incorporated  into  numerical 
models  and  presented  in  a  thesis  in  order  to  predict  the  characteristics 
of  aluminum  combustion  in  subscale  HTPB/Al/AP  propellant  motors. 
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Propeiiaat  Surface 


An  adequate  understanding  of  the  propellant  surface  phenojnena  ic 
probably  the  most  impoirtant  aspect  of  aluminum  combustion  modeling.  Its 
importance  is  due  to  the  fact  that  the  size  of  the  aluminum  and  aluminum 
oxide  particles  released  from  the  surface  largely  determines  the 
chemical  composition,  physical  properties,  and  the  distribution  of 
particles  throughout  the  rocket  motor  and  exhaust  plume.  The  complexity 
of  this  process  is  mainly  due  to  the  fact  that  the  aluminum  first  melts 
and  then  forms  ’'agglomerates''  that  can  become  .much  larger  than  the 
initial  particle  size  before  they  leave  the  propellant  surface  and 
ignite.  Both  the  size  and  degree  of  agglomeration  are  some  complicated 
function  of  the  ballistic  parameters  and  the  propellant  composition. 
There  have  been  several  attempts  to  model  this  phenomenon,  both 
experimental  and  analytical. 

Due  to  the  harsh  environment  in  the  combustion  chamber,  small 
aluminum  agglomerate  sizes (typically  <  100  microns),  and  relatively  fast 
combustion  times (typically  10  -  100  ms),  the  actual  agglomerate  sizes 
and  properties  are  difficult,  if  not  impossible,  to  observe  directly 
during  an  actual  motor  firing(Ref.  1).  Attempts  have  been  made  to 
measure  particle  sizes  by  collecting  the  quenched  exhaust  products  from 
burning  propellant  strands  in  closed  strana  combustion  bomb3(Ref.  2-4). 
By  collecting  the  exhaust  products,  scientists  were  able  to  measure  the 
pertinent  properties  of  the  combustion  products (i . e .  size,  composition, 
etc.).  Although  these  firings  should  have  produced  accurate  and 
consistent  measurements,  their  data  shows  significant  scatter.  These 
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experimental  methods  have,  however,  produced  some  very  accurate 
emperical  models  and  some  important  qualitative  inf ormation (Ref .  5). 

Another  experimental  method  that  has  been  used  to  determine  the 
combustion  processes  at  the  surface  of  the  propellant  is  the  use  of 
strand  burners  combined  with  high  speed  photography (Ref .  6-8) .  These 
usually  consist  of  burning  a  propellant  strand  in  a  pressurized 
combustion  bomb.  In  order  to  observe  the  surface  phenomenon,  the 
propellant  strand  is  sometimes  fed  into  the  combustion  bomb  at  a  rate 
equal  to  the  burning  rate  of  the  propellant.  In  this  way,  the 
propellant  surface  is  kept  stationary  relative  to  an  observation  window 
and  high  speed  film  is  used  to  record  the  combustion  processes  occurring 
at  the  propellant  surface.  The  use  of  these  combustion  bombs  has  also 
produced  some  qualitative  information  that  is  consistent  with  the 
exhaust  product  experiments.  It  should  be  noted,  however,  that  care 
should  be  taken  when  trying  to  extrapolate  Che  results  of  either  of 
these  experiments  to  actual  rocket  motors.  This  is  due  to  the  fact  that 
many  of  these  experiments  are  conducted  at  conditions  that  are  not 
consistent  with  actual  rocket  motor  environments (i.e.  low  chamber 
pressure,  different  species  concentrations,  etc.). 

As  noted  earlier,  both  of  these  experimental  methods  do  produce 
some  important  qualitative  information.  Some  of  these  dependencies 
include: 

-1)  as  chamber  pressure  increases,  the  agglomerate  size  decreases. 

-2)  as  the  normal  component  of  acceleration  increases,  agglomerate 
size  increases. 
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-3)  as  velocity  of  the  gas  increases,  agglomerate  size  decreases. 
-4)  as  aluminum  concentration  increases (or  oxidizer  concentration 


decreases),  agglomerate  size  increases. 

-5)  as  cunmonium  perchlorate (AP)  particle  size  decreases, 
agglomerate  size  increases. 

Another  method  that  can  be  used  to  model  the  propellant  surface 
phenomenon  is  to  apply  an  analytical  approach.  One  can  immediately  see 
the  difficulty  of  employing  an  analytical  model  just  by  loo)cing  at  the 
forces  that  occur  on  the  individual  agglomerates  as  they  emerge  from  the 
propellant  surface (see  Figure  1  and  Variable  Definitions (Ref .  9)).  The 
agglomerate  is  held  to  the  propellant  surface  by  surface  tension  and  a 
normal  acceleration  term(for  spin-stabalized  roclcets) .  The  forces  that 
attempt  to  "roll"  the  agglomerate  down  the  propellant  surface  include 
the  drag  and  axial  acceleration  while  the  lift  term  attempts  to  pull  the 
particle  away  from  the  propellant  surface.  Since  the  diameter  of  the 
agglomerate  is  a  function  of  time  and  the  flow  is  assumed  to  be 
turbulent,  ail  of  the  above  mentioned  forces  also  become  time  dependent. 
Add  to  this  the  fact  that  some  of  the  aluminum  may  be  converting  to 
aluminum  oxide  and  one  can  see  that  both  the  center  of  gravity  and 
center  of  pressure  of  the  particle  can  be  at  a  variable  location (causing 
additional  moments) .  It  should  be  noted  that  even  this  generalized 
force  diagram  has  made  several  simplifying  assumptions.  These 
assumptions  include:  spherical  agglomerate  shape,  smooth  propellant 
surface  and  axial  gas  velocity. 
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FIGURE  1:  GENERALIZED  FORCE  DIAGRAM 

(Ref.  9) 


Sn  - 
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VARIABLE  DEFINITIONS 


Ae  »  nozzle  exit  area 
A(x.t)  »  area  of  nozzle  or  plume  as  a 
function  of  axial  location 
and  time 

c.g.  =*  center  of  gravity 
c.p.  =  center  of  pressure 
D(t)  =  diameter  as  a  function  of  time 
E(x.r.t.T)  =  particle  emissivlty  as  a 
function  of  location,  time 
and  temperature 

Fd(x.r.t)  =  drag  force  as  a  function  of 
location  and  time 

FKx.r.t)  =  lift  force  as  a  function  of 
location  and  time 
gx(t}  =  axial  acceleration  as  a 
function  of  time 

gr(t)  =  radial  acceleration  as  a 
function  of  time 

m(t)  =  droplet  mass  as  a  function  of 
time 

R(x.t)  =  propellant  surface  location  as 
a  fwictlon  of  x  and  time 
Rt  =  radius  of  curvature  for  nozzle 
Sn  =  normal  surface  tension 
St  =  tangential  surface  tension 
T(x.r.t)  =  surface  shear  force  as  a 

function  of  location  and  time 
Ug(x.r.t)  =  gas  velocity  distribution  as 
a  function  of  location  and 
time 

Up(x.r.t.D)  =  particle  velocity  as  a 

function  of  location,  time 
and  diameter 
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One  type  of  analytical  approach  used  to  describe  the  propellant 
surface  phenomenon  involves  using  partial  differential  equations (Ref . 
10-11) .  First,  simplified  assumptions  about  the  propellant  surface 
processes  are  made.  Then,  the  necessary  conservation  equations  with 
their  appropriate  boundary  conditions  are  applied.  Finally,  the 
resulting  partial  differential  equations  are  usually  discretized  and 
incorporated  in  a  numerical  model.  The  accuracy  of  this  approach  is 
subject  to  the  initial  assumptions  made  and  the  step  size  taken  in  the 
discretization  scheme.  Depending  upon  the  simplifying  assumptions, 
these  methods  can  involve  considerable  computational  time  in  order  to 
produce  accurate  results. 

Another  analytical  method  employs  the  use  of  a  ’’pocket"  model  (Ref. 
12-14) .  Here,  the  propellant  is  assumed  to  be  made  up  of  a  series  of 
"pockets"  defined  by  increasing  AP  size  with  an  aluminum  particle  in  the 
center  of  the  pocket.  The  size  of  the  initial  aluminum  and  AP  determine 
if  the  aluminum  will  burn  immediately  or  agglomerate  and  then  burn  in 
the  next  larger  sized  pocket.  This  method  also  seems  to  produce 
accurate  qualitative  results  that  agree  with  experimental  data.  Since 
this  method  allows  for  varying  AP  size  and  can  be  easily  incorporated 
into  a  numerical  model,  it  is  probably  quicker  to  obtain  accuracy,  and 
therefore,  more  preferable  to  the  partial  differential  equation 
approach. 
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velocity  lags,  and  thermal  lags  for  that  matter,  they  usually  do  not 
consider  the  amount  of  slag  that  has  accumulated  or  the  corresponding 
mass  loss  to  the  flow (Ref  19).  This  may  be  due  to  the  fact  that  slag 
accumulation  is  largely  geometry  dependent.  For  this  reason,  only 
limited  studies  of  slag  accumulation  in  individual  solid  rocket  motors 
have  been  conducted (Ref .  20) . 

As  the  particles  enter  the  nozzle  throat  they  go  through  a  normal 
shock  wave.  This  increase  in  pressure  and  temperature  may  cause  many  of 
the  larger  particles  to  break  up  into  smaller  particles.  Some  models 
incorporate  this  fragmentation  by  considering  the  particle's  Weber 
nximber  and  comparing  it  to  a  critical  value  (Ref.  21).  Whenever  a 
particle  reaches  this  critical  value  the  particle  is  assumed  to  break 
up,  but  no  size  distribution  of  the  fragments  is  calculated.  In 
addition  to  the  fragmentation,  some  of  the  throat  material  may  ablate 
and  thus  add  "contaminants"  to  the  flow. 

As  the  particles  go  from  the  throat  to  the  exit  plane,  they 
encounter  numerous  oblique  shock  waves .  The  number  of  oblique  shock 
waves  is  usually  limited  by  careful  selection  of  the  nozzle  contour. 
These  oblique  shocks  may  also  act  to  shatter  the  larger  particles.  As 
the  particles  approach  the  exit,  some  of  the  smaller  particles  may  begin 
to  solidify.  Once  again,  the  smaller  particles  closely  follow  the 
streamlines  of  the  gas  while  the  larger  particles  lag  behind. 
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FIGURE  2:  GENERAL  NATURE  OF  ALUM  I  NUM  DROPLET 
COMBUSTION  IN  A  SOLID  ROCKET  MOTOR  ENVIRONMENT 
(Ref.  1) 
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is  recommendeci  that  further  investigation  of  the  combustion  processes  be 
conducted. 

Combustion  of  agglomerates  is  not  the  only  process  that  occurs  in 
the  combustion  chamber (see  Figure  3).  Fine  "smoke"  particles (usually 
submicron  in  size)  make  up  the  majority  of  the  combustion  products.  The 
agglomerates  may  break  up  due  to  internal  and  external  forces  or  they 
may  form  larger  agglomerates  due  to  collisions.  Both  of  these  competing 
processes  are  usually  considered  to  be  equal  and  no  net  change  in 
particle  sizes  is  calculated.  This  assumption  warrants  further 
investigation.  Also,  as  the  burning  particle  leaves  the  propellant 
surface,  it  emits  radiation  back  to  the  surface  and  to  other  cooler 
particles.  This  term,  too,  is  usually  neglected (even  though  it  has  been 
shown  to  contribute  as  much  as  twenty  percent  of  the  surface 
heating (Ref.  18)).  In  order  to  accurately  model  the  combustion  process, 
this  radiation  heat  transfer  must  be  included  in  both  the  surface  model 
and  the  combustion  chamber  model. 

Motor  Nozzle 

As  the  combustion  products  enter  the  nozzle  they  may  or  may  not  be 
completely  burned (see  Figure  4) .  As  the  agglomerates  enter  the  inlet 
the  smaller  particles (and  smoke)  will  follow  the  gas  streamlines  and 
velocities  more  closely.  This  lag  in  velocity  and  direction  by  the 
larger  particles  is  due  to  their  larger  mass  and  inertia.  For  this 
reason  and  due  to  the  geometry  of  the  nozzle,  some  of  the  larger  molten 
aluminum  and  aluminum  oxide  particles  may  accumulate  as  slag  on  the 
surface  of  the  nozzle  inlet.  Although  many  models  compensate  for  these 
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i  COMBUSTION  CHAMBER  PHENOMENON 


Particle  Agglomerate 
Formed  by  Collisions 


Ug(x.r.t) 


FIGURE  4:  NOZZLE  PHENOMENON 

“Smoke* 


Particles 


velocity  lags,  and  thermal  lags  for  that  matter,  they  usually  do  not 
consider  the  amount  of  slag  that  has  accumulated  or  the  corresponding 
mass  loss  to  the  flow (Ref  19) .  This  may  be  due  to  the  fact  that  slag 
accumulation  is  largely  geometry  dependent.  For  this  reason,  only 
limited  studies  of  slag  accumulation  in  individual  solid  rocicet  motors 
have  been  conducted (Ref .  20) . 

As  the  particles  enter  the  nozzle  throat  they  go  through  a  normal 
shocJc  wave.  This  increase  in  pressure  and  temperature  may  cause  many  of 
the  larger  particles  to  break  up  into  smaller  particles.  Some  models 
incorporate  this  fragmentation  by  considering  the  particle's  drag 
pressure  and  comparing  it  to  a  critical  value(Ref.  21).  Whenever  a 
particle  reaches  this  critical  value  the  particle  is  assumed  to  break 
up,  but  no  size  distribution  of  the  fragments  is  calculated.  In  this 
way,  however,  a  maximum  stable  droplet  size  can  be  calculated.  In 
addition  to  the  fragmentation,  some  of  the  throat  material  may  ablate 
and  thus  add  "contaminants"  to  the  flow. 

As  the  particles  go  from  the  throat  to  the  exit  plane,  they 
encounter  numerous  oblique  shock  waves.  The  number  of  oblique  shock 
waves  is  usually  limited  by  careful  selection  of  the  nozzle  contour. 
These  oblique  shocks  may  also  act  to  shatter  the  larger  particles.  As 
the  particles  approach  the  exit,  some  of  the  smaller  particles  may  begin 
to  solidify.  Once  again,  the  smaller  particles  closely  follow  the 
streamlines  of  the  gas  while  the  larger  particles  lag  behind. 
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If  one  assumes  a  uniform  particle  distribution  at  the  throat (which 
may  not  be  an  accurate  assumption),  as  in  most  models,  the  particle 
distribution  at  the  exit  consists  of  larger  particles  near  the  center 
and  smaller  particles  near  the  outside.  However,  the  nvmnber  of 
particles  along  the  centerline  at  the  exit  is  not  a  maximum.  The 
maximum  occurs  at  a  small  distance  from  the  centerline.  From  the 
maximum,  the  number  density  tapers  off  to  the  exit  radius  value (Ref. 

22)  . 

E.xhaust  Plume 

As  the  particles  enter  the  plume,  they  again  encounter  a  shoes 
wave  structure  consisting  of  Mach  disks (normal  shocks)  and  oblique 
shocks (see  Figure  5).  In  addition  to  these  shocks,  there  may  also  be  a 
bow  shock  at  the  edge  of  the  plume  depending  upon  the  velocity  and 
altitude  of  the  rocket.  The  plume  shape  itself  may  vary  drastically 
with  altitude - 

Another  phenomenon  which  occurs  in  the  plume  is  the  existence  of  a 
turbulent  mixing  layer.  This  mixing  layer  may  also  include  afterburning 
in  which  unburned  aluminum  particles  are  ignited.  Also,  some  of  the 
aluminum  and  aluminum  oxide  droplets  will  begin  to  solidify  in  the 
plume,  with  the  larger  droplets  solidifying  farther  downstream(Ref ,  23). 

A  large  portion  of  the  radiation  heat  transfer,  and  corresponding 
plume  signature,  comes  from  this  cloud  of  liquid  and  solid  aluminum 
oxide  particles.  Numerous  models  have  been  used  to  predict  the  plume 
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i  ROCKET  PLUME  PHENOMENON 


Particle  Radiation  to  Rocket  Body, 
Surroundings  and  Other  Particles 


radiation  using  Mie  theory,  Rayleigh  scattering,  and  other  methods (Ref. 
24-26) .  The  main  limit  to  the  use  of  these  radiation  models  is  an 
adequate  description  of  the  physical  properties  of  the  particles  and  an 
adequate  representation  of  the  size,  composition,  and  distribution  of 
the  particles  in  the  plume.  The  physical  properties  of  interest  include 
the  emissivity,  absorptivity,  reflectivity,  temperature,  velocity,  and 
density  of  the  particles.  Most  of  these  properties  are  inter-related 
and  some  of  them  have  been  measured  for  actual  rocket  motor 
environments (Ref .  27-28).  The  composition  of  the  particles  can  be 
greatly  affected  by  nozzle  ablation  and  afterburning.  In  turn,  the 
composition  of  the  particles  has  a  tremendous  influence  on  its  physical 
properties.  The  distribution  of  these  particles  is  mainly  a  function  of 
the  initial  agglomerate  size  and  the  nozzle  gecmetry. 

Goncluaians  and  ReconmendajLiaaa 

Although  some  very  accurate  numerical  models,  such  as  the  Solid 
Performance  Program(SPP  by  Hermsen  et.  al.)  and  the  One-Dimensional 
Reacting  Three-Phase  Flow  with  Mass  Transfer  Between  Phases (0D3P  by 
Kliegel  et.  al.),  exist,  it  is  recommended  that  some  improvements  should 
be  implemented (Ref .  29) .  These  improvements  include:  -1)  the 
incorporation  of  an  analytical  pocket  model,  which  includes  the 
influence  of  ballistic  parameters,  to  predict  the  size  and  extent  of 
surface  agglomeration. 

-2)  the  inclusion  of  radiation  heat  transfer  from  the  burning 
agglomerates  to  the  propellant  surface  and  other  particles. 
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-3)  an  adequate  model  to  compute  the  size  of  particles  formed  by 
collisions  of  agglomerates  in  the  combustion  chamber. 

-4)  a  description  which  includes  the  size  distribution  of  fragmented 
droplets  in  the  combustion  chamber,  the  nozzle,  and  the  exhaust  plume. 
-5)  a  simplified  relation  between  slag  accumulation  and  nozzle 
geometry. 

-6)  computation  of  the  actual  particle  distribution  at  the  nozzle 
throat . 

-7)  a  calculation  of  the  nozzle  ablation  rate  and  its  effect  on  plume 
radiation . 

Due  to  the  length  constraints  of  this  paper,  it  was  not  possible 
to  describe  any  area  of  the  rocket  motor  in  a  detailed,  analytical 
manner.  However,  some  of  these  models  and  their  suggested  improvements 
will  be  investigated  further.  Ideally,  these  improvements  will  be 
incorporated  into  a  numerical  model  and  presented  in  a  thesis  at  a  later 
date. 
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Abstract 

The  simulation  of  turbulent  and  Brownian  dispersion  of  solid  particles  in  an  atmospheric  boundary  layer 
requires  the  analysis  of  a  Lagrangian  trace  of  particle  trajectories.  A  computer  program  for  analyzing  the 
motion  of  solid  particles  in  the  turbulent  atmosphere  is  developed.  The  code  is  capable  of  providing  near¬ 
field  or  far-field  mass  concentrations  of  particles  from  continuous,  finite  duration,  and  instantaneous  point 
source  emissions.  The  fully  implicit  integration  of  the  particle  equation  of  motion  provides  particle  velocities 
induced  by  Stokes  drag,  Saffman  lift.  Brownian  diffusion,  and  gravity.  A  maximum  particle  concentration  of 
less  than  0.02%  by  weight  ensures  that  there  is  no  modification  of  the  air  flow  conditions  by  pjirticle  motion. 
Concentrations  of  this  order  allow  for  the  omission  of  all  particle-particle  interactions.  Three  sample  test 
cases  are  presented  for  illustrative  purposes 
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Introduction 


Problems  caused  by  particulate  dispersion  affect  our  lives  in  ways  we  hardly  even  realize.  The  eco¬ 
logical  difficulties  caused  by  acid  rain.  smog,  oil  spills,  and  forest  fires,  alter  the  environment  in  ways  that 
we  are  only  beginning  to  comprehend.  In  addition  to  these  and  numerous  other  natural  phenomena  involv¬ 
ing  particulate  transport,  there  are  many  man-made  tind  industrial  processes  for  which  particulate  flows 
play  an  essential  role.  For  instance,  in  the  early  1980’s.  the  microelectronics  industry  determined  that  the 
microcontamination  from  small  particles  is  the  leading  cause  of  loss  of  yield  in  the  manufacturing  process 
This  discovery  lead  to  a  massive  research  effort  in  the  area  of  microconlamination  control  and  clean  room 
applications.  In  addition,  it  has  been  determined  that  particles  deposited  from  the  air  flow  in  the  card 
passages,  limit  the  life  of  circuit  boards.  These  and  many  other  applications  lead  to  much  research  in  the 
area  of  aerosol  particle  dynamics. 

Due  to  the  earlier  research  efforts  in  the  area  of  aerosol  particle  dynamics,  a  wealth  of  knowledge  exists 
on  the  subject.  Fuchs  [1],  Hidy  and  Brock  [2],  van  de  Hulst  [3],  Twomey  (4],  and  Cadle  [5]  are  just  a  few 
of  the  classic  textbooks  on  the  subject.  Particle  dispersion  in  turbulent  flows  was  simulated  by  \hmadi 
[6],  Li  and  Ahmadi  [7],  and  Ounis  et.  al.  (8).  However,  to  the  author’s  knowledge,  no  detailed  Lagrangian 
simulations  of  atmospheric  turbulent  dispersion  have  been  published.  The  goal  of  this  study  is  to  provide  a 
large-scale  simulation  of  atmospheric  dispersion  of  fine  particles. 

The  use  of  a  supercomputer  facilitates  large  scale  simulations  which  may  provide  the  needed  insight 
into  the  processes  that  dominate  the  motion  and  dispersion  of  fine  particles  in  the  atmosphere.  This  report 
describes  a  simulation  program  currently  underway  at  Arnold  Engineering  Development  Center  for  simulating 
the  dispersion  of  fine  particles  in  the  atmosphere  with  application  to  environmental  compliance  with  the  Clean 
Air  Act  of  1990. 

Methodology 

In  this  simulation  a  steady-state  flow  field  is  taken  from  an  incompressible  Navier-Stokes  equation 
solver.  The  Lagrangian  equation  of  motion  for  a  heavy  particle  is  then  solved  in  time  for  a  release  of  particles. 
The  concentration  at  every  location  in  the  flow  is  taken  as  a  mass  percent  of  particles  in  a  given  region.  The 
distribution  of  particles  is  calculated  on  a  local  and  global  level.  Thirty-thousand  particles,  with  diameters 
ranging  from  0  to  100  microns,  are  simulated  for  instantaneous,  finite  duration,  and  continuous  releases. 

Mechanics  of  Aerosols 

By  definition,  an  aerosol  is  a  suspension  of  solid  or  liquid  particles  in  a  gas.  Common  aerosols  found 
in  nature  include:  dust,  smoke,  fog,  haze,  and  smog.  Aerosols  typically  range  in  size  from  0.001  to  100 
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microns  in  diameter.  To  place  these  dimensions  in  perspective,  note  that  the  mean  free  path  of  air  is  about 
0.07  microns  and  that  the  wavelength  of  visible  light  is  between  0.-4  to  0.7  microns  Also,  the  diameter  of 
0,01  microns  roughly  corresponds  to  the  transition  limit  between  molecule  and  particle.  Panicles  larger  than 
100  microns  typically  do  not  remain  suspended  in  air  for  a  significant  time  duration.  Thus,  the  behavior  of 
particles  is  significantly  affected  by  their  size. 


Aerosol  Particle  Motion 

The  Lagrangian  equation  of  motion  for  a  heavy  aerosol  particle  can  be  written  as 

duP  3irj^  j 
m—  =  -  «  )  +  ^9. 


(1) 


where  m  is  the  mass  of  the  particle  and  C«  is  the  Cunningham  correction  factor  given  by 


1.257  +  0.4 


exp 


Dividing  Eq.  (3.2)  by  infid/Ce  yields 

duf  .  f 
r-~  +  =  U-’  +  rp 

dt 

where  r  is  the  particle  relaxation  time.  This  is  defined  by 

_  mCe  _  d^p^Cc 
Zytfid  18/j 


where 


For  relatively  large  particles,  Cc  I  and 


m 


r  = 


dV 

18^ 


(2) 

(3) 

(4) 

(5) 

(6) 


Brownian  Motion 


Small  particles  suspended  in  a  fluid  undergo  random  translational  motion  due  to  molecular  collisions. 
This  phenomenon  is  referred  to  as  Brownian  motion.  The  Brownian  motion  leads  to  the  diffusion  of  particles 
in  accordance  with  Pick’s  law 


Ji  =  -D 


dc 

dxi 


(7) 


In  this  formula,  c  is  the  concentration,  J  is  the  flux,  and  D  is  the  diffusion  coefficient.  The  diffusivity  is 
given  as 

Jtr 

(8) 


kT 

D  =  zr^.Cc. 


Zirpd 

The  diffusion  coefficient  can  be  obtained  by  directly  substituting  Pick's  law  into  the  equation  for  mass 
conservation.  The  Brownian  forces  acting  on  the  particle  may  be  modeled  as  a  white  noise  process.  The 
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Brownian  motion  of  pariicies  is  modeled  by 


du 

—  +  du  =  n(f) 


19) 


with 

^  _  3ir^<f  _  1 

CfTTi  r ' 

Also,  n{<)  is  a  Gaussian  while  noise  process  with  a  spectral  intensity 


(10) 


^nn  —  ~  , 

Trm 

where  k  is  the  Boltzmann  constant  and  T  is  the  temperature.  From  Eq.  (9) 
spectrum  of  the  particle  velocity  is  given  as 


(11) 

it  follows  that  the  power 


=i  HiM)  |-  S„n(-) 


Here  the  system  function,  //(w),  is  defined  by 


/f(w)  =  - 


1 


Thus. 


S„u(w)  = 


iitj  *4"  0 

2kT0 


irm(w2  +  0'^) 


The  particle  autocorrelation  function  is  defined  as 


fiuu(r)  =  u'{t)u'{t  +  r), 

and  may  be  found  by  taking  the  inverse  Fourier  transform  of  the  power  spectrum,  i.e., 


The  diffusivity  is  then  given  by 


m 


Jo  Sir^id 


(12) 

(13) 

(14) 

(15) 

(16) 

(17) 


Algorithm  Development 

This  section  ts  concerned  with  the  actual  algorithm  used  in  the  LTIV13D  program.  The  lTIVIBD 
program  was  written  as  a  stand  alone  post  processor  for  the  Lawrence  Livermore  National  Laboratory 
FEMSA/B  gas  transport  and  dispersion  code.  The  LTIV13D  code  solves  the  equations  of  motion  for  soiio 
particles  undergoing  turbulent  and  Brownian  diffusion.  The  Stokes  drag  assumption  used,  limits  the  particle 
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Reynolds  number  to  a  maximum  of  1.  The  program  models  the  dispersion  of  randomly  sized  particles  using 
a  Lagrangian  approach.  The  computer  program  consists  of  a  mam  program  and  seven  subroutines  written 
in  standard  FORTRAN  “7  for  portability.  The  next  nine  sections  deal  with  the  specifics  of  the  code. 


Implicit  Integration 

The  main  program  performs  the  actual  integration  of  the  equations.  The  constants  that  are  used  and 
the  subroutines  that  are  called  are  described  in  this  section.  The  general  logic  is  as  follows.  The  program 
first  reads  the  nodal  points  and  velocities  from  an  external  file.  Typically,  these  files  are  output  from  the 
FEM3A/B  code.  The  new  particles  for  each  time  step  are  then  initialized.  For  each  particle,  the  mean  and 
the  fluctuation  fluid  velocities  are  found.  The  Brownian  force  is  also  evaluated.  Finally,  the  concentration 
is  calculated.  The  particle  velocity  and  position  are  then  updated.  Explicitly  the  particle  equations  in  Eq. 
(3.3)  are  given  as 

r— r—  -f  =  u-^  +  flu ,  ( 18) 

at 


T——  +  t/’=:v^-g+Bt.~q£.  (19) 

at 

4- +  flu).  (20) 

at 

Here  are  the  particle  velocities.  The  fluid  velocities  are  denoted  as  The  variables 

flu,  flu,  fltu  are  the  Brownian  forces,  E  is  the  electric  field  strength,  q  is  the  particle  charge  number,  t  is  the 
particle  relaxation  time,  and  g  is  the  gravitational  acceleration.  The  fully  implicit  formulation  uses  a  Crtmk 
-  Nicholson  discretization  method  given  by 


j.n  +  1  _ 


+  u" 


2  2  2 

A  Newton’s  iterative  method  is  used  for  the  solution  of  the  coupled  equations  in  u  and  r. 

The  same  routine  is  used  to  solve  for  these  quantities  for  the  next  particle,  trajecTocy  of  each  particle 
is  individually  advanced.  Because  of  the  amount  of  CPU  time  that  is  needed,  an  upoe-  :  n..i  of  30,000 
particles  is  suggested.  A  multiplicative  scale  factor  can  then  be  employed  to  simulate  a  t  ore  '■ea  stic  release 


Concentration  Calculation 

The  percent  by  mass  of  the  particle  in  the  area  surrounding  each  node  is  assigned  as  an  indication 
of  the  concentration  at  that  node  The  effect  of  mass  concentration  of  a  particle  on  each  mode  is 
inversely  proportional  to  the  distance  of  the  particle  from  that  node.  Specifically,  the  distance  of  the  nf^ 
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panicle  from  the  node  is  151  ven  as 


^nm  —  \/{^n  (Vn  l/m  )'  "*■  •^m)*  ^ 


where  dt  is  the  sum  of  the  distances  from  the  particle  to  its  8  surrounding  nodes,  is  the  mass  of  the 
n‘^  particle,  and  is  the  mass  of  the  fluid  in  the  element  containing  the  n‘^  particle.  The  mass  percent 
contribution  of  the  n‘*  particle  on  the  node  is  given  by 


^mn 


mP(d(  -dnm) 

m° 


(24) 


Turbulent  Fluctuations 

The  turbulent  velocity  field  is  constructed  by  the  superposition  of  a  fluctuating  velocity  onto  the 
steady  velocity  field.  A  Gaussian  random  variable  with  a  strength  u'  is  generated  in  each  of  the  three 
directions.  The  instantaneous  velocity  at  the  position  of  the  particle  is  given  as 

u  =  Zf+u'  (25) 

where  U  is  given  by  the  subroutine  fluid.  Then,  a  Gaussian  random  variable  with  strength  equal  to  the  rms 
velocity  may  be  used  as  a  fluctuating  component.  If  the  procedure  is  done  correctly,  an  imaginary  hot  wire 
held  at  a  specific  location  in  the  flow  will  see  a  signal  like  the  one  shown  in  figure  1.  from  Tennekes  and 
Lumley  [9]. 

Trilinear  Interpolation 

A  Newton's  iteration  method  and  a  trilinear  interpolation  (Benek  et  al.  [10])  is  used  to  find  the 
particle  position  in  grid  coordinates.  Then,  using  this  information,  it  interpolates  the  fluid  velocity  at  the 
particle  location.  Specifically,  suppose  the  nth  particle  located  at  the  point  xq,  is  contained  in  the  cell  with 
vertices 

Aj,  A*r,  A/  =  0, 1 

and  suppose  the  function  values 

/(*j+aj.t+a*,i+Ai).  Aj,  Afc,  Al  =  0, 1 

are  known.  The  value  /(xo)  is  then  approximated  as  follows.  First,  to  establish  a  relation  between  compu¬ 
tational  and  rectangular  coordinates  throughout  th**  cell,  define  the  trilinear  vector-valued  function: 

I 

X(q,/?,C)=  Y. 

F.1,f=0 

(26) 

=  «ooo  +  «ioo'}  +  <*010#^  +  «ooiC  +  crnoJ?/3 
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The  coefficients  ap,r  are  <ietermined  by  constraining  this  function  to  agree  with  known  coordinates  at  the 
cell  vertices: 


Q:p,r(Aj)^(  At)’(A/r  =  +  +  Aj,  Ai,  A/ -  0,  1.  (2i) 

p,},r  =  0 

This  is  a  linear  set  of  equations  which  can  be  solved  algebraically  for  the  components  of  the  vectors  ap^r- 
Similarly,  the  function  /{*)  is  represented  throughout  the  cell  with  a  trilinear  function: 

1 


F(r?,J,C)= 

p,q.r=0 


(28) 


=  /ooo  +  fiooV  +  /oiop?  +  /coiC  +  +  /lOlbC  +  /oil.^C  +  /un/^C- 

The  coefficients  fp^r  are  determined  by  constraining  this  function  to  agree  with  known  function  values  at 
the  cell  vertices; 

1 

Yi  /p»r(Aj)^(Alb)«(A/r  =  /(*j+A;,i+Ai,/+A/),  Aj,A*,A/  =  0,l.  (29) 

P.«,r=0 

Again,  this  is  a  linear  set  of  equations  which  can  be  solved  algebraically  for  /p,r.  To  find  the  computational 
coordinates  of  the  point  xo,  Newton’s  method  is  used  to  find  (r?o./^o.Co)  such  that  X{rio,  lio><^o)  =  *o- 
Specifically,  if  r?*  =s  (rik.litXk)  denote  successive  Newton  approximations  to  (r?o,/?o.Co).  then 

nt+i  =  >7t  -  {Df,X(T,t)]  [X(T;t)  -  xo]  (30) 


where 


n  Y  -  ^ 


dX  dX  dX 

dn'  dp  '  d<:  I 


(3i; 


After  this  iteration  converges  to  (tJo./^o.Co).  these  computational  coordinates  are  used  to  approximate  f[Xo] 
by  F(»7o,/?o,<o)- 


Results 

Each  test  case  uses  the  steady-state  flow  field  output  of  the  FEM3A/B  code.  Figure  2  shows  the 
three-dimensional  grid,  contours  of  the  average  u-component  of  velocity,  and  the  contours  of  the  average 
v-component  of  velocity  used  for  each  test  case.  The  45  x  12  x  15  grid  spans  a  three-dimensional  space 
of  300m  X  60m  x  120m.  The  grid  is  packed  near  the  ground,  and  in  the  area  of  the  release.  It  is  noted 
here  that  the  grid  used  in  these  cases  was  chosen  for  speed,  not  accuracy.  A  much  finer  grid  is  needed  for 
resolution  of  ail  relevant  scales.  The  acceleration  of  the  flow  by  the  ramp  is  a  classic  test  case.  The  results  of 
the  velocity  are  in  good  agreement  with  intuition  and  previous  data.  This  flow  was  chosen  because  intuition 
may  be  the  best  tool  available  for  solid  particulate  dispersion  concentrations  without  any  experimental  data 
for  comparison. 
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Instantaneous  Release 


In  the  first  test  case  30,000  panicles  are  released  simultaneously  from  a  point  source  in  the  flow.  Figures 
3  and  4  depict  the  diffusion  dynamics  of  the  instantaneous  release.  Figure  3  shows  the  initial  distribution 
of  the  particle  diameters,  the  contours  of  initial  concentration,  and  the  spatial  distribution  of  the  particles. 
The  linear  distribution  between  0  and  100  microns  has  an  average  of  50  microns  and  a  standzurd  deviation 
of  25  microns.  Approximately  300  particles  of  a  given  size  exist  in  the  flow  field.  The  initial  concentration 
is  centered  at  the  point  source.  The  concentration  is  highest  at  the  source  and  decreases  monotonically  in 
a  sphere  surrounding  the  source.  The  average  diameter  of  50  microns  reflects  that  the  particles  are  not 
separated  with  size,  but  instead  are  well  mixed. 

Figure  4  shows  the  dust  cloud  87.5  seconds  after  the  release.  The  distribution  is  now  truncated  at  50 
microns.  This  indicates  that  the  larger  particles  have  precipitated  out  of  the  flow.  The  mean  diameter  is 
now  23  microns  and  the  standard  deviation  is  17  microns.  The  small  particles  are  still  suspended  in  the 
flow.  The  concentration  87.5  seconds  after  release  shows  that  the  mass  percent  is  centered  at  the  middle  of 
the  cloud  because  of  the  relatively  large  number  of  small  particles  as  compared  with  the  smaller  number  of 
large  particles.  The  cloud  has  moved  downstream,  leaving  only  deposited  particles  on  the  ground.  Almost 
no  concentration  is  found  in  the  wake  of  the  instantaneous  released  cloud.  The  cloud  has  diffused  somewhat 
in  the  vertical  direction  due  to  Brownian  and  turbulent  fluctuations.  The  spatial  distribution  of  the  particles 
left  in  the  flow  show  that  the  smaller  particles  are  much  further  downstream  than  the  larger  particles.  The 
larger  particles  have  a  greater  relaxation  time  and  take  longer  to  respond  to  the  fluctuations  of  the  velocity. 
The  larger  particles  also  occupy  the  space  only  near  the  j"ound.  Oi-nvity  acts  strongly  on  these  particles 
and  forces  precipitation.  The  smaller  particles  extend  quite  further  from  the  ground  because  of  the  reduced 
effect  of  gravity. 

Finite  Duration  Release 

In  the  second  test  case,  400  particles  are  released  every  half-second  for  37.5  seconds  from  a  point  source 
in  the  flow.  This  gives  a  total  release  of  30,000  particles  analogous  to  test  case  1.  Figures  5  and  6  depict 
the  diffusion  dynamics  of  the  finite  duration  release.  Figure  5  shows  the  initial  distribution  of  the  particle 
diameters,  the  initial  concentration  contours,  and  the  spatial  distribution  of  particles  at  time  tg.  The  linear 
distribution  between  0  and  100  microns  has  an  average  of  50  microns  and  a  standard  deviation  of  25  microns. 
Approximately  15  particles  of  a  given  size  exist  in  the  flow  field  initially.  The  relative  scatter  of  the  diameters 
has  now  increased  because  of  the  lower  number  of  samples.  The  initial  concentration  is  centered  at  the  point 
source.  The  concentration  is  highest  at  the  source  and  decreases  monotonically  in  a  sphere  surrounding 
the  source.  For  the  initial  condition,  the  particles  occupy  the  region  surrounding  the  source.  The  average 
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diameter  of  50  microns  reHecis  that  the  particles  are  not  separated  with  size,  but  instead  are  homoeeneously 
mixed. 

Figure  6  shows  the  dust  cloud  87.5  seconds  after  the  release.  This  figure  indicates  that  the  distribution  is 
now  constant  for  the  first  50  microns  then  decreeises  to  80  microns.  The  largest  particle  left  in  the  flow  field 
is  shown  to  be  80  microns.  This  indicates  that  the  larger  particles  have  precipitated  out  of  the  flow  The 
smaller  particles,  however,  are  still  suspended  in  the  flow  even  from  the  initial  release.  The  concentration 
variation  shows  the  cloud  extending  far  downstream  of  the  initial  release.  The  mass  percent  is  centered 
near  the  trailing  edge  of  the  cloud  because  of  the  relatively  large  number  of  large  particles  still  in  the  flow 
as  compared  with  the  instantaneous  release.  The  cloud  of  the  finite  duration  release  is  much  longer  than 
the  cloud  of  the  instantaneous  release.  The  cloud  has  moved  downstream  leaving  only  deposited  particles 
on  the  ground.  Almost  no  concentration  is  found  in  the  wake  of  the  finite  duration  cloud.  The  cloud  has 
diffused  somewhat  in  the  vertical  direction  due  to  Brownian  and  turbulent  fluctuations.  Here,  we  also  see 
that  the  smaller  particles  are  much  further  downstream  than  the  larger  particles.  The  smaller  particles 
released  initially  have  now  traveled  far  downstream  and  the  larger  particles  released  at  the  end  of  the  finite 
duration  have  not  yet  precipitated  out. 

Contur  rovis  Release 

In  the  third  test  case,  150  particles  are  released  every  half-second  for  100  seconds,  from  a  point  source 
in  the  flow.  This  gives  a  total  release  of  30,000  particles  which  is  analogous  to  test  cases  1  and  2.  Figures 
7  and  8  depict  the  diffusion  dynamics  of  the  continuous  release.  Figure  7  shows  the  initial  distribution  of 
the  particle  diameters,  the  initial  concentration,  and  the  spatial  distribution  of  the  particles  at  time  to-  The 
linear  distribution  between  0  and  100  microns  has  an  average  of  50  microns  and  a  standard  deviation  of  25 
microns.  Approximately  5  particles  of  a  given  size  exist  in  the  flow  field  initially.  The  relative  scatter  of 
the  diameters  has  again  increased  because  of  the  number  of  samples.  The  initial  concentration  is  centered 
at  the  point  source.  The  concentration  is  highest  at  the  source  and  decreases  monctonically  in  a  sphere 
surrounding  the  source.  For  the  initial  condition,  the  particles  occupy  the  region  surrounding  the  source. 
The  average  diameter  of  50  microns  reflects  that  the  particles  are  not  separated  with  size,  but  instead  are 
homogeneously  mixed.  This  is  again  very  similar  to  both  the  instantaneous  and  finite  duration  releases. 

Figure  8  shows  the  dust  cloud  87.5  seconds  after  the  initial  release.  This  figure  indicates  that  the 
distribution  is  now  constant  for  the  first  50  microns  then  decreases  to  100  microns.  The  largest  particle 
left  in  the  flow  field  is  shown  to  be  100  microns.  The  large  particles  still  exist  in  the  flow  because  they 
are  released  with  every  time  step.  The  relative  number  of  large  particles  indicates  that  the  larger  particles 
have  precipitated  out  of  the  flow.  The  small  particles,  however,  are  still  suspended  in  the  flow  even  from 
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I  he  muial  release.  I'lie  ('oncentration  shows  the  eloud  extendiiii;  far  <iownsireani  of  the  uutial  release  I  h< 
mass  percent  is  centered  near  the  trailing  edge  of  the  cloud  because  of  the  relatively  large  number  of  large 
particles  still  in  the  flow.  The  continuous  release  cloud  is  now  very  long  because  of  the  particles  in  the  flow 
field  from  every  time  step  The  concentration  decreases  rapidly  with  downstream  distance.  Here,  we  see  that 
the  smaller  particles  are  much  further  downstream  than  the  larger  particles.  The  small  particles  released 
initially  have  now  traveled  far  downstream  and  the  larger  particles  released  at  the  last  time  step  have  not 
yet  precipitated  out.  The  monoionic  decrease  in  the  average  particle  diameter  is  almost  logarithmic  with 
downstream  distance. 

Conclusions 

The  main  conclusion  drawn  here  is  the  availability  of  a  usable  tool  for  the  simulation  of  solid  particulate 
diffusion  and  transport  in  atmospheric  turbulent  flow,  instantaneous,  finite  duration,  and  continuously 
released  particles  have  been  simulated  for  randomly  dispersed  particles.  This  indicates  great  hope  for  the 
eventual  pollution  predictive  capability. 

Future  Work 

A  logical  improvement  for  future  capability  is  to  integrate  the  LTM3D  program  and  the  FEM3A/B 
program  to  solve  for  gaseous,  liquid,  and  solid  particulate  dispersion  simultaneously.  The  inclusion  of  tem¬ 
perature  and  density  variations  could  also  be  integrated.  The  combustion  and  chemical  reaction  of  solid 
particles  also  may  need  to  be  accounted  for. 
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Figure  1:  Gaussian  Distributed  Random  Variable 
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Figure  2:  Flow  Geometry  and  Velocities  for  Test  Cases 
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Concent rjt ;on 


Figure  3;  Initial  Condition  for  Instantaneous  Release 
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Figure  8:  Cloud  Orientation  87.5  seconds  After  Initial  Release 
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CAO  AND  ACOUSTIC  BEM  APPLIED  TO  THE  MODELLING  OF  THE  AEDC  A8TF  EGM8 
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Abeliact 

CAD  input  and  display  packages  were  interfaced  to  acoustic  boundary  element  codes.  These  codes  were  examined 
with  respect  to  accuracy  of  amplitude,  phase,  and  frequency,  as  well  as  the  rates  of  convergences  as  functions 
of  dement  resolution.  Sound  pressure  levd  distributions  were  collated  at  two  resolution  levels  for  geometries 
and  boundary  conditions  corresponding  to  certain  AEDC  ASTF  EGMS  diffusers.  A  fulhscale  acoustic  boundary 
element  modd  of  the  AEDC  ASTF  segment  containing  the  EGMS  was  constructed. 


4-2 


CAD  AND  ACOUSTIC  BEM  APPLIED  TO  THE  MODELLING  OF  THE  AEOC  A8TF  EGMS 


Rkhard  A.  Maradwll 


INTRODUCTION 

in  1991,  under  the  guidance  of  R.  R  Jonee  III  of  the  Sverdrup  Technology  lnc./AEDC  Group,  M.  A.  Weaver  [1] 
began  to  inveatigate  the  feasibility  of  using  a  threo-dimensional  acoustic  boundary  element  code  to  predict  the 
acoustic  response  of  the  Arnold  Engineering  Development  Center's  (AEDC)  Aeropropubion  Systems  Test  Faciiity 
(ASTF)  Exhaust  Gas  Management  System  (EGMS).  A  description  of  the  AEDC  ASTF  EGMS  can  be  found  in 
the  report  by  Weaver  [1]  as  vvefl  as  a  history  of  the  problem  which  will  not  be  repeated  here.  Weaver  [1]  reported 
that  for  geometries  similar  to  the  EGMS  diffuser  section,  an  acoustic  boundary  element  method  (BEM)  using 
isoparametric  elements  seemed  to  predict  resonant  frequencies  and  modal  amplitude  distributions,  at  least  along 
the  boundaries.  This  report  describes  a  continuation  and  extension  of  the  previous  study.  In  particular,  this  study, 
also  under  the  direction  of  R.  R.  Jones  III,  attempted  with  some  degree  of  success  to  scale  up  the  use  of  an 
acoustic  boundary  element  method  to  realistic  full-size  problems. 

OBJECTIVES  ACHIEVED 

1.  Interfaced  a  three-dimensional  computer-aided  drawing  package  to  the  construction  of  the  acoustic  boundary 
element  code  (BEMAP)  input  data  files.  Now  both  the  acoustic  boundary  element  and  field  point  inputs 
can  be  graphically  prepared  and  displayad. 

2.  Interfaced  a  three-dimensional  high  resolution  color  CAD  data  display  package  to  the  outputs  of  the  acoustic 
boundary  element  code  BEMAP. 

3.  Examined  boundary  element  acoustic  cedes  HELM2D1  (developed  at  the  Univereity  of  Tennessee)  and 
BEMAP  (developed  at  the  Univeraity  of  Kentucky)  with  respect  to  accuracies  of  amplitude,  phase,  and 
frequency,  as  weH  as  the  rates  of  convergence  as  functiorw  of  element  resolution. 

4.  Examined  BEMAP  behavior  at  two  resolution  levels  for  geometries  and  bouttdary  conditions  with  known 
analytical  resulta.  In  particular,  computed  sound  pressure  level  (SPL)  distributions  over  a  150  Hz  frequency 
range,  were  compared  to  theoretical  and  experimental  results. 

5.  Comtructed  a  full-scale  acoustic  boundary  element  model  of  the  AEDC  ASTF  segment  containmg  the  EGMS 
with  a  cylindrical  diffuser  section.  This  model  contained  196S  linear  acoustic  boundary  elements  and  over 
1700  field  pointa. 

OBJECTIVE  NOT  ACHIEVED 

Difficulties  with  the  BEMAP  code  and/or  its  interaction  with  the  computational  platform  on  larger  problems 
prevented  the  computation  of  a  full-scale  model  over  the  frequency  bands  of  interest  during  the  short  duration  of 
this  contract. 

FIGURES  AND  DISCUSSION 

Figures  1-21  foHow  on  page  4,  discussion  begira  on  page  19. 
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Figure  1:  toerior  lighted  view  of  coane  acoustic  boundary  element  discretization  of  a  cylindrical  model  of  the  Arnold 
Engineering  Development  Center  (AEOC)  Aeropropulsion  Systems  Test  Facility  (ASTF)  Exhaust  Gas  Man¬ 
agement  System  (EGMS).  The  axis  is  along  the  i-axis.  The  center  of  the  left  endeap  is  at  the  origin. 
Vwous  boundary  conditions  were  applied  at  the  endcaps.  the  circumference  of  the  cylinder  was  taken  as 
perfectly  ngid.  Acoustic  sources  were  modelled  both  within  the  interior  and  on  the  left  endcap  for  various 
confi^rations.  The  coarse  models  used  198  linear  acoustic  boundary  elements. 
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Cutaway  interior  tii^ted  view  o(  the  coarse  acoustic  boundary  element  discretization  of  the  cylindrical  model 
of  the  AEDC  ASTF  EGMS.  For  the  cases  where  an  offset  endpanel  source  was  used,  it  was  modelled  bv 
specifying  a  velocity  on  the  “middle'  panel  of  the  fi:st  quadrant  {^-z  plane)  of  the  left  endcap.  This  view 
particularly  shows  how  rough  an  approximation  a  12-sided  polygon  is  to  a  cylindrical  surface.  Not  surprisingly, 
these  models  had  difficulty  reproducing  circumferential  modes. 


Figure  3:  Extefiof  lii^hted  vww  of  ihe  fine  >co<ntic  boundary  element  diecretization  of  the  cylindrical  model  of  the 
AEDC  A5TF  EGMS.  Note  that  compared  to  the  coarae  model  displayed  in  Figures  1  and  2,  the  role  of  the  r 
and  r  axes  are  reversed.  Here,  the  r  axis  is  in  the  axial  direction.  This  was  to  facilitate  boundary  condition 
editing  in  the  fine  model  setups.  These  models  consisted  of  1104  linear  acoustic  boundary  elements.  Using 
the  BEMAP  code,  each  frequency  need  about  10  hours  to  compute  on  a  1-2  MFLOP  computer. 
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Figure  4:  Surface  view  of  'he  plar>es  of  interiof  field  node  points  used  in  both  the  coarse  and  fine  cylindrical  acousiic 
models  of  the  AEDC  ASTF  EGMS.  The  axes  orientation  shown  here  corresponds  to  the  coarse  boundarv 
element  models.  The  fine  model  field  points  are  the  same,  except  the  x  and  r  axes  are  reversed.  The  interior 
field  point  nodes  arc  at  the  intersections  of  the  grid  lines.  The  number  of  field  points  computed  was  the 
same  lor  both  the  coarse  and  fine  models  (195  points).  The  field  point  "endcap  planes  were  inset 
of  the  evlinder  length  from  the  boundary  clement  model  endcaps.  The  inset  distance  in  the  radial  direction 
was  similar. 
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Filjure  6:  Coarse  model,  offset  panel  source.  =  -  •  H 


Closest  to  exact  resonant  /  =  IT  -  Hz. 


Shows  start  ol  rrsonattce  somewhat  hieh 
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Firure  9:  Coant  model,  offset  p#nel  source.  ^  Hz 
Peak  of  computed  resonance. 


Slightly  off  peak.  2.8  Hz  away  from  exact. 
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Sli|?ht  indication  o<  neai  resonance. 
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Figure  13; 


Figure  14: 
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Figure  16:  Fine  model,  offset  panel  source,  '  =  !  i: 
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Well  o»T  resonance. 
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f  igure  17a-17b  :  The  order  of  linear  acoustic  boundary  element  methods  increases  with  frequency,  although  the  error  itself 

becomes  larger.  That  is,  relative  to  the  dimensions  of  the  problem,  at  low  frequencies  increasing  the  number 
of  elements  improves  on  an  already  good  solution  only  slightly,  at  mid-frequencies  the  improvement  becomes 
quite  significant  (about  2nd  order),  and  at  high  frequencies  the  improvement  is  dramatic  (4th  order  or  higher). 


amplitude  errors  versus  number  op  elements 

(Scaling  13  arbtlrory,  votues  deoeno  on  frequency, 
oroblem  dimensions,  element  sires,  and  other  toctors 
Results  are  typical  for  relatively  high  frequencies. 

1^-  i.e.  about  three  boundary  elements  per  wavelength.) 
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Figure  13:  Exterior  lighted  view  of  the  CHAMBERI  acoustic  boundary  element  model  o)  the  AEDC  ASTF  segment 
containing  the  EGMS.  This  view  shows  the  downstream  location  at  the  top  of  the  page.  i.e.  ilow  is  in 
the  +e-direction.  The  exit  endcap  shown  has  an  acoustic  impedance  boundary  condition  corresponding  to 
perfectly  absorbing  to  normally  incident  waves.  The  two  circumferential  rings  of  elements  nearest  the  exit 
endcap  also  have  this  boundary  condition.  This  model  consists  of  1963  linear  acoustic  boundary  elements. 
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Figure  19;  Cutaway  lighted  view  q(  the  CHAMBERl  acoustic  boundary  element  model  o<  the  AEDC  ,  ,STF  segment 
containing  the  ECM5.  The  ddTuser  region  opening  into  the  mam  tunnel,  and  the  bulkhead  region  are  ciearlv 
shown.  The  interior  cylindrical  diffuser  region  has  the  same  dimensions  and  element  resolution  as  the  fine 
resolution  cylinder  model. 
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Figure  20:  Cutaway  lighted  y-z  view  of  the  CHAMBERl  acoustic  boundary  element  model  of  the  AEDC  ASTF  segment 
containing  the  EGMS.  The  cylindrical  diffuser  has  a  relatively  high  boundary  element  resolution  in  order  to 
capture  interior  resonance  frequencies  and  modal  amplitude  distributions  accurately.  The  coarser  downstream 
tunnel  segment  discretization  is  believed  to  be  adequate  for  modelling  the  diffuser  exit  acoustic  impedance. 
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Pigure  31;  Field  point  planes  for  the  CHAMBERl  acoustic  boundary  element  model  depicted  in  Figures  15-20.  Field 
points  are  at  the  intersections  of  grid  lines.  As  in  the  cylindrical  models,  the  field  point  planes  are  inset 
a  small  distence  from  the  boundary  element  surfaces.  The  field  point  resolution  here  is  higher  than  the 
cylindrical  models  with  over  1748  field  points. 
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CAD  lnttai»ce 

Although  perhaps  more  accurate  than  linear  acoustic  boundary  elements,  the  isoparametric  acoustic  boundary 
elements  of  the  earlier  study  [1]  have  the  rather  severe  disadvantage  of  beii^  almost  impossible  to  interface  to 
CAD  draering  and  display  packages.  In  order  to  use  acoustic  boundary  element  methods  as  a  practical  desipi 
tool  for  modelling  realistic  geometries,  such  CAD  interfaces  are  a  necessity.  Linear  acoustic  boundary  elements 
are  not  only  much  easier  to  interface,  the  models  are  also  easier  to  debug  since  all  nodes  we  at  the  vertices  of 
the  elements,  i.e.  all  nodes  are  at  the  grid  line  intersectkms.  For  these  reasotw,  this  study  used  linear  acoustic 
boundary  elements  exclusively. 

Aoougtk  Boundmy  Ekmatt  L<m  RttMolution  Modda 

An  acoustic  boundary  element  model  of  a  cylindrical  diffuser  section  of  the  AEDC  ASTF  EGMS  with  a  minimal 
number  of  elements  was  constructed  to  test  BEMAP  with  linear  quadrilateral  elements  and  to  compare  to  the 
earlier  results  by  Weaver  [1].  The  lew  resolution,  or  coarse  model,  had  about  three  elements  per  wavelength  in  all 
directkms  at  the  highest  frequencies  of  interest.  The  coarse  model  geometry  is  shown  in  Figures  1  and  2.  Note 
that  the  coarse  model  here  has  about  the  same  number  of  elements  as  the  *1100  element  model  of  the  earlier 
study  [1].  The  fine  or  high  resolution  model  of  this  study,  Figuie  3,  would  have  been  very  difficult  to  enter  into 
BEMAP  without  the  use  of  a  CAD  input  interface. 

The  source  location  and  boundary  conditiofis  (hardwaii  everywhere)  were  repeated  and  similar  results  were 
obtained,  although  here  the  resonant  peaks  were  less  well  defined.  That  is,  the  isoparametric  elements  used  by 
Weaver  [1]  produced  sharper  resorances  somewhat  closer  to  the  exact  analytical  frequencies.  The  present  study 
solved  for  interior  field  points,  whereas  the  previous  study  examined  the  boundaries  only.  The  field  point  locations 
are  depicted  in  Figure  4. 

In  addition  to  the  comparisons  with  earlier  study  results,  computer  runs  were  made  with: 
s  Weaver's  source  location,  but  zero  pressure  endcap  boundary  conditions, 
s  Source  point  at  the  center  with  zero  pressure  endcap  conditions. 

s  Offset  wall  panel  source  with  the  remainder  of  that  endcap  having  a  zero  velocity  or 
hardwaU  condition,  and  the  other  endcap  having  a  zero  pressure  boundary  condition. 

A  comparison  of  this  last  set  of  conditions  at  low  and  high  resolutions  is  presented  in  this  report. 

Acoustic  Boundmy  Element  High  Reaohtion  Modah 

The  high  resolution  or  fine  acoustic  boundary  element  geometry  used  to  model  the  cylindrical  diffuser  portion  of 
the  AEDC  ASTF  EGMS  is  shown  as  Figure  3.  These  modek  with  1104  linear  quadrilateral  acoustic  boundary 
elements  and  the  fidd  point  geometries  of  Figure  4  were  run  for  a  variety  of  boundary  conditions.  Principal  among 
these  boundary  conditions  were: 

s  Source  point  in  the  center  with  zero  pressure  endcaps. 

s  Source  point  in  the  center  with  one  endcap  vebcity  (hardwaii)  and  the  other  zero  pressure. 

•  Offset  wall  panel  source  with  the  remainder  of  that  eiKlcap  having  a  zero  velocity  (hardwaii) 
condition,  and  the  other  endcap  having  a  zero  pressure  boundary  condition. 
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Comtpmmon  of  the  Two  Reoohitiom 

A*  M  iSuttesbv*  eumpia  of  the  charactar  of  retuftt  obtainod  at  the  two  reaoiutions,  the  aound  preaaurc  lewd 
(SPL|  diatributiom  were  commuted  arouitd  a  known  analytical  rcaonance  at  85.2  Hz  for  the  case  of  the  offaet  wall 
pand  tource  and  boundary  conditiorw  of  one  hardwwall  endcap ,  the  other  zero  preeaure. 

Fijuraa  5-10  ahow  the  SPL  diatribution  in  dB  for  the  reaolution  modd  r  ound  the  exact  85.2  Hz  reaonance. 
The  cartaapendinf  raauHa  for  the  high  reaolution  modd  are  diapfayed  aa  Fip.jrea  11-15.  Briefly  comparing  the  two 
aeta  of  laaulta  ahow  that  the  high  reaolution  modd  haa  ita  reaorvance  doaar  to  the  exact  value,  and  that  it  ia  much 
aharper.  The  low  reaolution  modd  dkplaya  the  reaonance,  but  at  a  higher  in  frequency  than  predicted  analytically. 
Theae  raauita  are  typkd  when  the  low  reaolution  medd  doea  in  fact  capture  the  reaonance.  Often,  however,  the 
IcMr  reaolution  modd  miaaee  reaonanoea  entirely.  Thia  effect  wee  noted  by  Weaver  (1]  for  Moparametric  etementa 
end  ia  true  for  linaar  quadriiataral  elcmenta.  It  appeara  that  aix  or  more  linear  acouatic  bouitdary  dementa  per 
wavelangth  in  aU  diroctiona  are  needed  to  rdiably  capture  reaonancea. 

Further  EWeeta  of  Efement  Reeohitiott 

Exammation  of  the  literature  revealed  almoet  no  publiahed  reaulta  on  the  element  reeolutiona  required  to  aolve 
varioua  actud  acouatic  problema.  Reference  [3]  recommended  at  leaat  four  dementa  per  wavelength,  and  adviaed 
increaaing  the  number  of  dementa  until  the  “errof  waa  aufficiently  amali.  An  extenahre  inveatigation  of  the  effecta 
of  dement  reeohitien  on  the  amplitude  end  phaee  errora  in  bouitdary  elemont  methoda  waa  undartakan  uaing 
the  code  HELM201  for  the  two>dimenamnal  Helmholtz  equation.  Trenda  obaerved  in  two-dimeneioiia  ware  then 
examined  in  three  danenaiona  uaing  BEMAP  to  tea  if  they  held  true  in  general.  Except  for  an  abaolute  acale  factor 
which  me  dependent  on  the  varioua  other  aapecta  of  the  code  implenientationa  and  problem  particuian,  a  number 
d  traita  laemed  to  generalize.  Among  the  moat  interaadng  were  thoae  ahown  in  Figurca  17a  and  17b,  that  ia, 
the  order  of  liiteer  acouatic  boundary  dement  methoda  increaaca  with  frequency,  although  the  error  itaetf  becomea 
larger.  In  other  worda,  relative  to  the  dimenaione  of  the  problem,  at  low  fiequendea  increaaing  the  number  of 
dementa  improvea  on  an  already  good  aolution  only  alightiy,  at  mid-frequendea  the  improvement  becomea  quite 
a^nificant  (about  2nd  order),  end  at  high  ftequendea  the  improvement  ia  dramatic  (4th  order  or  higher). 


Acouatic  boundary  dement  methoda  continue  to  hold  promiae  to  be  able  to  naodd  the  interior  sound  diatrdMition 
and  reaonancea  of  the  AEDC  ASTF  EGMS.  A  code  written  apecificaliy  for  this  purpose  could  take  advantage 
of  pertacular  characteristica  of  the  geometry,  frequeiKy  ranges,  and  accuracy  requirements  in  oriler  to  arrive  at 
reaulto  within  realistic  wall-clock  times. 
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A  MULTIGRAPH  IMPLEMENTATION  OF  A 
DISTRIBUTED  IMAGE  PROCESSING  SYSTEM 


Michael  S.  Moore 
PhD.  Candidate 

Department  of  Electrical  Engineering 
Vamderbilt  University 


Abstract 


At  Arnold  Engineering  Development  Center,  videos  of  rocket  plumes  are  used  in  anal¬ 
ysis.  The  videos  are  inherently  noisy,  so  they  must  be  di^tally  processed  before  they  are 
useful.  Processing  the  video  frames  on  normal  digital  computers  can  require  days  or  even 
weeks.  Thus,  there  is  a  need  for  a  high  speed  image  processing  system.  It  has  been  seen 
that  the  speed  of  image  processing  operations  can  be  greatly  increased  by  distributing  the 
computational  load  over  several  workstations,  PCs,  or  transputers  using  Multigraph.  Multi¬ 
graph,  a  system  integration  tool  developed  at  Vanderbilt  University,  allows  the  building  of 
complex  algorithms  from  simpler  processing  blocks.  Muitigraph  is  capable  of  distributing 
processes  over  a  network  to  a  variety  of  computer  architectures,  and  the  network  does  not 
have  to  be  homogeneous. 

This  report  presents  the  results  of  the  1992  Graduate  Summer  Research  Program:  a 
Multigraph  implementation  of  a  distributed  image  processing  system.  The  system  utilizes  a 
non- homogeneous  network  of  workstations  to  gain  a  considerable  increase  in  the  execution 
speed  of  image  enhancement,  noise  reduction,  and  analysis  algorithms.  The  flexibility  and 
speed  of  the  system  have  been  demonstrated.  The  overall  results  of  the  summer  research 
are  very  promising. 
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introduction- 


The  engineers  at  Arnold  Engineering  Development  Center  (AEDC)  routinely  use  rocket 
test  videos  to  analyze  rocket  plumes.  The  videos  are  noisy,  so  digital  image  processing 
techniques  are  used  to  reduce  noise  and  enhance  the  digitized  video  sequences.  Due  to  the 
huge  amount  of  data  involved  '  and  the  complexity  of  the  algorithms,  it  can  take  days  or 
even  weeks  to  process  the  video  from  one  test  firing.  Thus,  there  is  a  need  for  a  low  cost, 
flexible,  high  speed  image  processing  system  with  which  the  engineers  can  reduce  greatly 
the  computer  time  and  man-hours  needed  to  process  the  video  data. 

Another  need  is  for  a  system  with  which  engineers  can  rapidly  develop  new  image 
processing  algorithms.  It  is  desirable  that  the  user  be  able  to  build  an  algorithm  from 
a  library  of  small  pre-coded  operation  blocks  without  knowing  the  internal  details.  The 
system  should  be  based  on  models  that  allow  the  algorithms  to  be  easily  modified  without 
recompiling  of  the  source  code. 

This  work  represents  a  continuation  of  a  research  effort  initiated  by  this  author  as  a  par¬ 
ticipant  of  the  1991  AFOSR-RDL  summer  research  program  [4],  The  new  distributed  image 
processing  system  has  been  greatly  improved  and  new  functions  have  ’’een  added.  The  sys¬ 
tem  is  designed  to  speed  up  image  processing  routines  by  distributing  the  computational 
load  across  a  network  of  Unix  workstations.  The  machines  used  include  Sun  SPARC2,  IBM 
RS6000,  and  4d340  Silicon  Graphics  Iris  stations.  Since  all  of  these  machines  are  commonly 
used  at  AEDC,  the  distributed  image  processing  system  requires  no  new  hardware.  The 
Multigraph  execution  environment  is  used  for  data  flow,  process  control,  and  communica¬ 
tion.  The  HDL  (hierarchically  descriptive  language)  interface  is  used  to  define  graphs  of 
each  of  the  processing  algorithms  and  to  communicate  with  the  Multigraph  kernel. 

In  the  following  text,  an  explanation  of  the  objectives  of  the  summer  research  is  given. 
Then  a  brief  description  of  the  tools  used  in  implementing  the  distributed  image  processing 
system  is  discussed.  The  image  processing  capabilities  of  the  system  are  discussed  next,  with 
a  brief  description  of  each  of  the  utility  functions  and  processing  algorithms  that  have  been 
implemented  to  date.  Then  the  performance  of  the  system  is  analyzed  and  recommendations 
for  future  research  are  made. 

’For  example  consider  a  standard  60  second  video;  60  Irames/sec  interlaced  (net  30  frames/sec),  512  by 
480  pixels,  one  byte/pixel.  That  is  more  than  442  Mbytes  of  data  to  be  processed  per  minute  of  video! 
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OBJECTIVES  OF  THE  RESEARCH  EFFORT: 


The  basic  research  objectives  were: 

•  Create  image  processing  utilities  to  perform  tasks  such  as  reading,  writing,  and  passing 
images.  Create  a  basic  image  viewer  for  on  line  use. 

•  Establish  protocol  for  system  topology  and  image  structure  to  be  used  in  future  ad¬ 
ditions  to  the  library  of  image  processing  functions. 

•  Begin  a  library  of  general  image  processing  routines.  Create  an  interface  between  the 
new  protocol  and  previously  written  Multigraph  routines.  [6] 

•  Improve  upon  image  splitting  and  reconstruction  techniques  so  as  to  reduce  the  effects 
of  a  non-homogeneous  network. 

•  Demonstrate  the  speed  and  flexibility  of  the  system  with  several  examples  of  dis¬ 
tributed  implementations  of  image  processing  algorithms. 

As  wUl  be  described  in  the  text,  these  objectives  were  met  with  the  summer’s  work. 


THE  MULTIGRAPH  ENVIRONMENT: 

Multigraph  is  a  system  integration  tool  developed  at  Vanderbilt  University.  It  allows 
the  construction  of  complex  algorithms  from  smaller  modules  of  code.  Muitigraph  will  au¬ 
tomatically  schedule,  execute,  and  control  the  data  flow  of  an  algorithm.  Since  Multigraph 
supports  distributed  processes,  any  part  of  an  algorithm  can  be  executed  remotely.  The 
remote  processes  can  be  executed  on  workstations,  PCs,  or  transputers,  and  the  network 
need  not  be  homogeneous. 

The  heart  of  Multigraph  is  the  Multigraph  kernel  (MGK).  It  does  the  scheduling,  process 
control,  and  the  communication.  It  uses  a  model  of  the  system  called  a  graph.  A  graph 
represents  all  aspects  of  the  algorithm  to  be  implemented.  Important  components  of  the 
graph  are  given  below. 


5-4 


1.  ^cforsare  the  rompulationaJ  operators  of  the  graph.  To  each  actor  is  attached  a  state, 
which  can  be  inactive,  active,  ready,  or  running.  The  state  is  used  by  the  scheduler 
to  decide  when  the  actor’s  function  will  be  performed  [1]. 

An  actor  is  made  up  of  several  components: 

(a)  The  script  is  the  actual  code  that  ’will  be  executed  when  the  actor  runs.  The 
actor  script  can  be  in  almost  any  language,  and  can  be  used  by  any  number  of 
unique  actors  in  the  same  graph.  Within  the  script.  Multigraph  kernel  calls  are 
used  to  receive  the  input  data,  propagate  the  output  data,  and  manage  memory. 

(b)  The  context  of  the  actor  is  an  area  of  static  memory  allocated  by  the  kernel  that 
is  not  reset  or  cleared  each  time  the  actor  is  run.  Thus,  the  actor  can  save  its 
state  between  runs.  Also,  the  context  may  be  set  when  the  graph  is  built  to  pass 
initial  parameters  to  the  script.  This  type  of  context  allows  the  execution  to  be 
controlled  by  the  specification  of  the  graph  with  no  need  to  recompile  the  script. 

(c)  The  input  ports  are  where  data  is  passed  to  the  actor  by  the  kernel  and  the  output 
ports  are  where  data  leaves  the  actor.  The  number  of  input  and  output  ports  of 
an  actor  is  specified  when  the  actor  is  created  and  cannot  be  changed  without 
rebuilding. 

(d)  The  control  principle  is  the  method  by  which  the  kernel  decides  when  to  schedule 
an  actor  for  execution.  An  actor  can  be  an  ifall  or  an  ifany  actor.  An  ifall  actor 
will  be  scheduled  for  execution  when  all  of  its  input  ports  contain  data.  An 
ifany  actor  will  be  scheduled  whenever  one  of  more  of  its  input  ports  have  data 
available. 

2.  Datanodes  are  places  in  memory  where  data  is  queued  between  actornodes.  They 
connect  actor  output  ports  to  actor  input  ports.  Since  they  are  queued,  they  allow 
the  kernel  to  provide  dynamic  scheduling  for  the  graph.  A  datanode  can  be  connected 
to  any  number  of  actor  input  ports,  and  any  number  of  actor  output  ports  can  be 
connected  to  it.  (Note  that  data  cannot  flow  directly  from  an  actor’s  output  port  to 
another’s  input  port.  It  must  flow  through  a  datanode.  This  is  because  the  input 
and  output  ports  do  not  queue  data.)  Each  datanode  is  assigned  a  type  -  either 
scalar  or  stream.  When  an  actor  receives  data  from  an  input  port  attached  to  a  scalar 
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datauode,  only  the  last  piece  of  data  propagated  to  the  datanode  will  be  received. 
However,  data  on  a  stream  datanode  will  be  received  in  FIFO  order. 

3.  The  Environment  to  which  an  actor  is  attached  is  a  protected  set  of  system  resources 
that  it  uses.  An  environment  has  associated  with  it  a  priority,  which  the  kernel  uses 
in  scheduling.  Only  one  actor  per  environment  can  be  executed  at  any  time.  Thus, 
for  all  actors  associated  with  a  particular  environment,  mutual  exclusion  is  ensured 
on  that  environment’s  resources  . 

4.  A  task  is  an  interface  to  the  physical  resources  of  the  machine.  In  a  multitasking 
environment  tasks  are  the  different  processing  threads  available  to  the  Multigraph 
kernel.  In  a  multiprocessor  system  they  are  the  individual  processors  themselves  [1]. 
The  system  discussed  in  this  work  is  a  multiprocessor  system,  so  the  tasks  represent 
the  workstations  that  will  actually  do  the  computation.  Each  environment  is  attached 
to  a  task,  and  any  task  can  have  multiple  environments.  Also,  multiple  tasks  can  be 
associated  with  a  machine. 

The  Multigraph  kernel  provides  the  interface  to  build,  modify,  monitor,  and  control 
the  graph.  ^  When  constructing  large  graphs,  though,  the  direct  MGK  interface  calls  can 
become  cumbersome.  To  allow  more  elegant  communication  with  the  kernel,  higher  level 
tools  exists.  One  such  tool  is  the  HDL  (hierarchal  descriptive  language)  interface.  The  HDL 
interface  is  described  in  the  following  section. 


THE  HDL  INTERFACE: 

HDL  (hierarchal  descriptive  language)  is  used  to  model  the  signal  flow  graph  of  the 
system.  The  language  represents  the  processing  system  in  terms  of  its  structure.  The  HDL 
language  is  naturally  analogous  to  a  block  diagram  representation  of  the  graph,  and  as  such 
provides  an  easy  method  of  describing  the  processing  blocks  and  interconnections  of  complex 
systems.  As  its  name  infers,  HDL  is  inherently  hierarchal.  Thus,  any  block  of  a  graph  can 
represent  an  arbitrary  hierarchy  of  sub-blocks.  This  allows  for  high  level  abstraction  and 
simpler  data  flow  graphs.  The  simplest  and  lowest  level  blocks  are  called  primitives.  They 
*For  a  complete  description  of  the  MGK  interface,  see  [2]. 
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are  the  blocks  to  which  actual  code  is  associated.  Blocks  made  up  of  one  or  more  sub-blocks 
are  called  compounds.  The  sub-blocks  of  a  compound  can  be  either  primitives,  compounds, 
or  any  mixture  of  each. 


The  Semantics  of  HDL 

The  HDL  language  uses  scripts  with  a  specific  syntax  to  interpret  the  structure  of  the 
system.  Either  a  primitive  or  a  compound  is  declared  in  each  script.  Examples  of  a  primitive 
and  a  compound  script  are  ^ven  in  figures  1  and  2  respectively.  Note  that  HDL  ignores 
anything  after  a  on  a  line,  so  lines  starting  with  are  comments.  Referring  to  figure 
1,  note  that  the  primitive  is  named  mclean  primitive.  It  is  an  IFALL  actor  with  two  inputs 
and  one  output.  Its  context  consists  of  a  static  parameter  called  mcleanjnumber,  which  is 
an  integer  that  defaults  to  0  at  build  time  if  it  is  not  otherwise  set.  The  context  could  also 
have  dynamic  parameters,  but  does  not.  None  of  the  actors  in  this  image  processing  system 
use  dynamic  parameters.  See  [3]  for  a  complete  description  of  dynamic  parameters.  This 
actor  is  aittached  to  the  environment  named  “e”,  which  is  attached  to  the  task  named  “t”. 
The  script  that  will  be  executed  when  “mclean  primitive”  is  run  is  named  “mclean”.  The 
script  “mclean”  must  have  been  linked  to  the  Multigraph  kernel  at  compile  time. 

The  compound  definition  shown  in  figure  2  contains  the  primitives  OpenRas  primitive^ 
SetParams  primitive,  mclean  primitive,  and  WriteRas  primitive.  It  has  no  inputs,  outputs, 
static  puameters,  or  dynamic  parameters.  The  datanodes  of  a  graph  are  specified  by 
the  SIGNALS  declaration  in  the  script.  The  data  nodes  of  the  mclean  compound  are 
mcleanjparameters,  topology,  rawjmage,  and  new  image.  The  meanings  of  the  PARAMS, 
SHARED,  VARS,  and  COMPUTE  declarations  can  be  found  in  the  HDL  manual  [3],  but 
are  not  important  here.  The  STRUCT  declaration  is  where  the  topology  of  the  actor  is 
defined.  Each  sub-actor  has  input  and  output  ports  which  are  attached  to  data  nodes.  Each 
sub-actor  is  also  attached  to  an  environment  and  a  task.  Note  that  in  figure  2  no  static 
parameter  is  passed  to  the  mclean  primitive  (ie.  NIL  is  in  the  static  parameter  position). 
In  this  situation,  the  graph  builder  will  assign  the  default  value,  which  is  0,  to  the  mclean 
pnnutive's  local  static  parameter  mclean..number. 

To  build  the  mclean  graph,  the  HDL  scripts  for  the  mclean  compound  and  all  of  the 
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Definition  for  loclean  primitive 
(DEFPRIMITIVE  loclean  primitive!  :IFALL 
;  input  and  output  definitions 
( 

( loclean.parametersl  rSTREAH)  (iRasInl  .-STREAM)  ->  (iRasOut!  rSTREAM) 

) 

;  static  peurameters 
(  (  Imclean.number I  0  :INT  )  ) 

;  dynamic  parameters 
NIL 

; environment /task  mapping 
(lei)  (jtl) 

;  t«o  unused  graphics  slots 
NIL  NIL 
; script  name 
"mclean" 

) 


Figure  1:  Example  of  an  HDL  Primitive 

related  primitives  must  be  loaded  using  the  HDL  interface.  For  specific  procedures  for 
loading  scripts,  building  graphs,  and  executing  algorithms,  see  [3]. 


TOPOLOGICAL  CQNSIDEP.,ATTQNS: 

For  any  distributed  system  special  considerations  must  be  addressed.  A  computational 
topology  must  be  chosen  for  the  algorithm.  In  other  words,  one  must  decide  how  the 
processing  is  to  be  broken  up  between  the  available  resources.  The  best  choice  is  highly 
dependent  upon  both  the  algorithm  and  the  resources  available. 

One  way  to  distribute  the  work  load  is  to  break  the  algorithm  into  steps,  or  computa¬ 
tional  blocks.  Each  computational  block  can  be  assigned  to  a  processor,  thus  presumably 
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;;;  hdl  compound  ’mclean’ 

(DEFCOHPOUND  I mclean I  <  ->  )  NIL  NIL 

;no  inputs,  outputs,  static  or  d3mamic  parameters 
; environment /task  mapping 
del)  (Itl) 

;  teo  unused  graphics  slots 
NIL  NIL 

(SIGNALS  (Imclean.parametersl  :STREAM)  (Itopologyj  -.STREAM) 

(IraB.imagel  ;STREAM)  (Ines.imagel  rSTREAM)  ) 

(PARAMS)  (SHARED)  (VARS)  (COMPUTE) 

(STRUCT  (I  Open  I  (lOpenRas  primitive! 

(  ->  Iraw.imagel 
NIL  NIL 

del)  (Itl)  )  ) 

dSetParamsI  (iSetParams  primitive  I 

(  ->  Imclean.parametersl  I  topology  I 
NIL  MIL 

del)  dtl)  )  ) 

(I morphological  cleaner]  ( I mclean  primitive  I 

(  Imclean.parametersl  I  raw.  image  I  ->  ] new. image  I ) 
NIL  NIL 

del)  dtl)  )  ) 

( I  Write  I  (iVriteRas  primitive  I 
(iRasOutj  ->  ) 

NIL  NIL 

del)  dtl)  )  ) 

) 

) 


Figure  2:  Example  of  an  HDL  Compound 
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Input  data 


output  data 


Figure  3:  A  Graph  of  a  Decomposed  Algorithm 

causing  a  net  speedup  in  processing.  See  figure  3.  However,  note  that  if  one  processing 
block’s  input  data  depends  upon  another  block’s  output  data,  as  in  part  a  of  figure  4,  then 
the  system  is  elfectively  a  serial  processor.  Processor  n  will  have  to  wait  on  processor  n  —  1 
to  finish  before  it  can  run,  Tliis  serial  topology  can  still  be  effective  if  the  data  set  is  such 
that  many  “waves”  of  data  wiU  be  propagated  through  the  system  one  after  the  other. 
Then,  if  each  of  the  serial  processing  blocks  is  assumed  to  have  approximately  the  same 
execution  time,  then  each  processor  “down  the  line”  will  be  idle  only  until  the  first  wave  of 
data  reaches  it.  After  that,  a  continuous  flow  of  data  will  be  achieved.  Note  that  the  pro¬ 
cessing  blocks  are  assumed  to  have  similar  execution  times.  The  difficulties  of  this  condition 
become  clear  in  light  of  the  fact  that  the  processor  assigned  to  block  n  could  much  faster 
or  slower  than  the  one  assigned  to  block  n  -  1.  Since  most  networks  are  non- homogeneous, 
this  is  a  likely  situation.  Also,  an  algorithm’s  natural  decomposition  may  not  lead  to  blocks 
of  similar  complexity,  which  would  similarly  effect  the  system’s  speed. 

Some  algorithms  may  decompose  naturally  in  a  parallel  manner,  as  show  in  part  b 
of  figure  3.  Again,  the  successful  speedup  depends  upon  the  condition  that  the  parallel 
processing  blocks  are  executed  at  nearly  the  same  speed.  Thus  all  of  the  inputs 

to  block  m  arrive  with  approximate  synchronicity.  This  will  insure  that  block  m  does  not 
waste  time  waiting  on  one  or  more  input(s)  while  its  others  are  already  available.  Again, 
non-homogeneity  of  the  processor  network  becomes  an  issue. 
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processor  N-1  — <{  processor  N 


Part  a:  Serial  Processing  Topology 


Part  b:  Parallel  Processing  Topology 


Figure  4:  Serial  and  Parallel  Decomposition 

A  different  approach  to  parallelizing  a  system  is  to  distribute  the  data  instead  of  dis¬ 
tributing  the  algorithm.  In  the  image  processing  case,  this  means  to  split  the  image  into 
pieces  and  send  the  pieces  to  remote  processors.  Each  remote  processor  will  perform  the 
same  operation  on  its  sub-image.  The  processed  sub-images  are  then  put  together  to  form 
the  output  image.  This  technique  is  commonly  referred  to  as  farming.  Not  all  imaging 
algorithms  are  appropriate  for  farming  because  they  do  not  use  neighborhood  operations. 
Neighborhood  operations  form  the  output  at  each  pixel  from  pixels  in  a  neighborhood 
around  that  position.  Examples  of  neighborhood  operators  are  the  median  filter,  the  mor¬ 
phological  opening  and  closing  operations,  and  linear  convolution  with  a  short  duration  FIR 
filter.  By  choosing  sub-images  with  slightly  overlapping  edges,  a  distributed  implementation 
of  a  neighborhood  operation  can  have  the  same  outcome  as  the  non-distributed  case.  (ie. 
Neighborhood  operations  can  be  considered  invariant  under  this  type  of  farming  process.) 

Farming  was  the  process  distribution  method  used  in  implementing  all  of  the  algorithms 
presently  supported  by  the  distributed  imaging  system.  However,  the  other  methods  can 
easily  be  applied  where  appropriate.  An  example  of  an  applicable  algorithm  is  the  2- 
dimensional  wavelet  transformation.  The  planned  approach  to  parallelizing  the  Daubechies 
wavelet  transform  for  images  includes  a  combination  of  farming  and  process  decomposition 
methods. 
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IMAGE  PROCESSING  UTILITIES: 


Several  utility  actors  were  developed  that  are  used  by  almost  all  of  the  implemented 
algorithms.  Basic  I/O,  system  topology  control,  and  image  viewing  capabilities  were  needed. 
The  following  Multigraph  actors  were  written  to  fit  these  needs. 

1.  OpenRas  is  the  actor  that  opens  an  image  file,  reads  it  into  Multigrapb  memory,  and 
propagates  it  to  an  output  port.  The  type  of  memory  structure  it  allocates  is  called 
Raslmg.  The  Rasimg  structure  is  described  in  [4]  and  is  declared  in  the  file  include.h. 
OpenRas  is  always  capable  of  reading  compressed  or  uncompressed  Sun  rasterfiles, 
whether  or  not  compression  is  specified  in  the  filename.  If  PBMplus  is  available  on 
the  computer  where  OpenRas  is  executed,  it  can  read  any  PBM  supported  image 
format,  compressed  or  uncompressed.  In  its  present  state,  OpenRas  runs  once  and 
then  deactivates  itself.  This  property  can  be  eaisily  modified  to  accommodate  the 
processing  of  image  batches. 

2.  WriteRas  receives  an  image  buffer  and  puts  it  into  a  Raslmg  type  structure  if  it 
recognizes  the  inconaing  buffer  type  It  then  writes  an  uncompress  Sun  rasterfile. 

3.  SpliiRas  receives  an  image  and  a  two  data  structures  containing  information  describing 
the  topology  of  the  process  distribution  and  the  size  of  the  neighborhood  that  will  be 
used.  From  this  information  it  decides  how  to  cut  the  image  and  how  much  overlap 
to  allow  along  the  cuts.  It  then  propagates  the  sub-images  from  its  outputs. 

4.  ConstructRas  receives  the  topology  information,  the  width  and  height  of  the  original 
image,  and  the  sub-images.  It  reconstructs  an  output  image  of  the  same  size  as  the 
original  image,  throwing  away  the  overlapping  edges.  The  result  is  propagated  to  an 
output  port.  An  interesting  feature  of  ConstructRas  is  that  it  is  an  IFANY  actor. 
Each  time  one  of  its  inputs  becomes  available,  its  script  is  executed.  The  output 
image  buffer,  which  is  stored  as  context,  is  updated  whenever  a  sub-image  becomes 
available.  The  partially  filled  out  image  is  propagated  from  an  output  port  called 

^  A  m&gic  number  specifies  each  supported  buffer  type.  Each  actor  in  the  system  checks  the  magic  number 
of  the  image  buffer  it  receives  and  tries  to  convert  the  image  to  its  format.  These  conversions  are  built  into 
the  image  receiving  functions,  so  they  need  not  be  apparent  to  the  script  writer 
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display.update  to  the  image  displayer.  Thus,  the  user  can  evaluate  the  results  as  the 
output  image  is  still  being  formed. 

5.  SetTopology  allows  control  of  the  image  splitter.  The  operator  can  choose  the  number 
of  slices  to  be  made  in  the  horizontal  and  vertical  directions.  At  the  present  time,  the 
maximum  total  number  of  sub-images  is  10.  This  constraint  is  set  only  by  the  HDL 
scripts,  not  by  the  actor  scripts.  It  can  easily  be  made  as  large  as  desired  by  editing 
the  HDL  code. 

6.  DisplayUpdate  receives  an  image  and  displays  it  in  an  X- window.  The  displayer  was 
adapted  from  xloadimage,  which  is  a  standard  X-windows  image  viewer. 

Along  the  way,  an  image  structure  protocol  was  adopted  so  that  future  efforts  will  be 
compatible  with  the  present  system.  Interfaces  between  already  existing  structures  and 
the  new  protocol  were  developed.  The  supported  formats  axe  IMAGE  Rasimg,  and 
AEDC_EL3  ®.  Any  actor  that  uses  one  of  these  image  data  structures  internally  can  receive 
a  generic  image  by  calling  its  interface  function.  (For  example  Receive  Rasimg  receives  an 
image  and  puts  it  into  a  Rasimg  type  buffer.) 


IMAGE  PROCESSING  LIBRARY: 

The  following  is  a  brief  description  of  the  algorithms  presently  available  on  the  dis¬ 
tributed  image  processing  system.  These  algorithms  were  chosen  to  be  implemented  first 
because  they  are  computationally  expensive  and  thus  usually  require  large  amounts  of  pro¬ 
cessing  time.  Many  other  algorithms  can  easily  be  added  in  the  future. 

1.  Mclean  is  a  morphological  noise  reduction  technique  developed  by  R.  A.  Peters  of 
Vanderbilt  University  [5].  It  is  based  on  the  morphological  operations  open  and  close. 
Mclean  is  extremely  effective  in  removing  noise  from  images  while  leaving  edges  and 
small  grained  features  intact.  This  algorithm  was  implemented  during  the  1991  RDL 
summer  research  program  by  this  author  [4].  Slight  changes  have  been  made  to  the 

’The  IMAGE  data  structure  is  defined  in  image.h 

^The  AEDC_EL3  was  adopted  by  the  EL3  group  at  AEDC  as  the  standard  Multigraph  image  format  to 
be  used  in  the  future.  It  is  defined  in  elSJmage.h 
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software  since  then,  but  the  algorithm  ib  still  as  described  in  {4]  and  [5],  The  mciean 
actor  receives  a  data  structure  called  mclean.paramelerJist,  which  is  created  by  an¬ 
other  aolor  called  SeiPammeiers.  The  mclean_parameter  Jist  contains  the  structuring 
element  specifications  and  other  information  that  mciean  needs  (see  [4]).  The  mciean 
parameters  can  be  set  to  default  or  edited  at  run-time,  so  the  algorithm  can  be  ad¬ 
justed  without  recompiling  the  code. 

2.  The  adaptive  histogram  equalization  (AHE)  algorithm  is  usuadly  used  for  edge  and 
detail  enhancement.  The  algorithm  rescales  the  intensity  map  for  each  pixel  based  on 
the  histogram  of  a  neighborhood  around  it.  Note  that  the  effect  of  the  AHE  operation 
is  highly  dependent  upon  the  neighborhood  size.  In  general,  large  neighborhoods  will 
tend  to  bring  out  details  such  as  edges  very  well.  The  actor  was  written  by  Rich 
Souder,  and  uses  the  IMAGE  data  structure  [6].  Only  minor  adjustments  were  made 
to  the  code,  including  the  installation  of  the  data  structure  conversion  interface  pre¬ 
viously  discussed.  The  AHE  actor  receives  width  and  height  of  the  neighborhood  via 
a  data  structure  called  neighborhood Jnfo,  which  is  created  by  the  actor  SetTopology. 

3.  The  median  filter  script,  which  w'as  also  written  by  Rich  Souder,  implements  a  neigh¬ 
borhood  based  image  smoothing  technique.  The  algorithm  makes  a  list  of  intensities 
that  occur  in  a  window  around  each  pixel.  It  sorts  the  list,  and  sets  the  pixel  in  the 
center  of  the  window  equal  to  the  median  value.  The  effect  is  similar  to  that  of  an 
averaging,  or  rolling  ball  filter.  The  noise  is  reduced  in  proportion  to  the  window 
area.  However,  unlike  the  mciean  filter,  the  edges  and  small  grained  features  are  sig¬ 
nificantly  degraded.  The  median  filter  algorithm  receives  the  same  neighborhoodJnfo 
data  structure  that  the  AHE  actor  does,  so  the  window  size  is  easily  changed. 


SYSTEM  PERFORMANCE: 

The  parallel  image  processing  system  was  tested  extensively.  As  an  example,  the  results 
of  the  parallel  mciean  system  test  will  be  given  here.  Note  that  the  other  algorithms 
implemented  have  graphs  very  similar  to  the  mciean  graph  shown  in  figure  5,  and  they 
performed  similarly  in  tests. 


5-14 


A  graphical  model  of  the  parallel  mclean  system  is  shown  in  figure  5.  It  uses  a  compound 
called  lOjmcleans,  which  contains  10  unique  mclean  actors.  Each  of  these  mclean  actors 
could  be  assigned  to  a  separate  processor,  but  they  do  not  have  to  be.  Two  of  more  of  them 
can  be  executed  on  the  same  processor  if  desired. 

The  mclean  system  was  developed  and  tested  on  a  network  of  seven  Unix  workstations. 
A  diagram  of  the  network  is  given  in  figure  6.  ®  Note  that  each  station  has  attached  to  it 
at  least  one  task,  and  that  each  task  has  an  associated  environment.  The  master  machine, 
named  nemo,  has  two  tasks-environment  pairs.  The  master  environment  performs  the  image 
I/O,  splitting,  and  reconstruction.  The  io_u)in  environment  sets  the  parameters  and  runs 
the  displayer.  The  slave  environments  are  attached  to  the  mclean  actors  and  therefore  do 
most  of  the  computation. 

The  speedup  data  shown  in  figure  7  w’ere  taken  during  the  1991  Summer  Research 
Program.  See  [4].  The  system  was  implemented  on  a  network  of  Sun  SPARCl  workstations. 
The  algorithm  tested  was  a  parallel  version  of  the  mclean  filter  which  used  from  one  to 
three  remote  processes.  Note  that  the  execution  time  of  the  mclean  algorithm  is  highly 
dependent  upon  the  structuring  element  size.  Thus,  the  structuring  element  size  (actually 
structuring  element  width)  was  varied  as  well  as  the  number  of  remote  processes.  Each 
test  was  performed  on  a  512  by  512,  8-bit  greyscale  image.  Note  that  when  the  structuring 
element  size  varies,  the  execution  time  seems  to  increase  in  a  quadratic  fashion.  This  is 
because  of  the  number  of  operations  involved  in  mclean  is  proportional  to  the  area  of  the 
structuring  element,  which  is  the  square  of  the  parameter  we  are  calling  structuring  element 
size. 

The  speedup  due  to  the  paraUelization  was  considerable.  As  an  example,  refer  to  the 
second  plot  of  figure  7.  The  speedup  figures  for  the  curve  marked  X  =  23  are  given  by 


Sp€edup[2]  = 
Speedupid]  = 


465seconds 

26bseconds 

46Sseconds 


=  1.75 


=  2.27 


(1) 


(2) 


205seconds 

where  Spe€dup[n]  is  the  speedup  ratio  observed  with  n  remote  processes.  Note  that  the 
speedup  ratio  is  less  than  n  due  to  the  overhead  associated  with  data  communications  and 


'Other  configurations  were  tested  that  used  from  two  to  ten  workstations.  The  diagram  shows  the  most 
commonly  used  configuration. 
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MCLEAN  SYSTEM  SCHEMATIC  (version  2) 
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Figure  5:  Graphical  Model  of  the  Parallel  mclean  System 
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the  non-homogeneity  of  the  network.  The  speedup  ratios  are  expected  to  improve  with 
a  future  modification  of  the  splitting  and  reconstruction  actors.  These  awrtors  will  soon 
be  improved  to  adaptively  rout  data  based  on  the  apparent  performance  of  the  processing 
nodes.  This  will  reduce  the  effects  of  using  a  non-homogeneous  network. 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH: 

In  the  future  much  work  should  be  done  to  improve  the  Multigraph  distributed  image 
processing  system.  The  recommended  improvements  and  areas  of  study  include* . 

•  Additions  to  the  library  of  image  processing  routines  to  include  a  broad  range  of 
common  applications. 

•  Improvement  of  the  image  data  distribution  technique  to  include  adaptive  routing  of 
data. 

•  A  general  study  of  topology  and  algorithm  parameterization.  This  study  should  con¬ 
sider  the  problem  of  matching  processing  topologies  to  a  classes  of  image  processing 
problems.  The  goal  should  be  to  develop  an  analytic  method  matching  the  optimal 
system  configuration  to  a  given  algorithm. 
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The  rounded  off  rectangular  boxes  connected  to  the  machines  represent  tasks. 
The  rectangular  boxes  represent  environments.  Each  machine  is  a  either  a  Sun 
SPARC,  an  IBM  RS600,  or  a  4d340  Silicon  Graphics  workstation. 


Figure  6:  System  Development  and  Test  Network 


5-18 


EXECUTION  TIME  IN  SECONDS  EXECUnON  TIME  IN  SECONDS 


500 


EXECUTION  TIME  vi.  S.E  SIZE 


References 


[1]  Abbott,  B.A.,  Bapty,  T.A.,Biegl,  C.:  “Experiences  Using  Model-Based  Techniques  For 
The  Development  Of  A  Large  Parallel  Instrumentation  System”,  Final  Report  for  the 
1992  USAF-RDL  summer  research  program,  1992. 

[2]  Biegl,  Csaba.:  “Multigraph  Kernel  (MGK)  User’s  Manual”,  Dept,  of  Electrical  Engi¬ 
neering,  Vanderbilt  University,  1988. 

[3]  Karsai,  G.:  “Hierarchical  Description  Language  (HDL)  User’s  Manual”  Dept,  of  Elec¬ 
trical  Engineering,  Vanderbilt  University,  Technical  Report  #87-004,  1987. 

[4]  Moore,  Michael  S.:  “Multigraph  Implementation  of  Image  Morphology”,  Final  report 
for  the  1991  USAF-RDL  summer  research  program,  1991, 

[5]  Peters,  Richard  Alan.:  “Image  Sequence  Noise  Reduction  using  Morphological  Filters”, 
Final  report  for  the  AFOSR.  Research  Initiation  Program,  1991. 

[6]  Souder,  Richard  S.:  “Parallel  Distributed  Image  Processing”,  Master’s  Thesis  Submit¬ 
ted  to  the  Faculty  of  the  Graduate  School  of  Vanderbilt  University,  1989, 


5-20 


A  CELL  AVERAGED  APPROACH  TO  THE  SOLUTION 
OF  INTEGRAL  CONSERVATION  LAWS 


Blair  H.  Rollin 

Graduate  Research  Assistant 
Department  of  Mathematics 


The  University  of  Tennessee 
Space  Institute 
Tullahoma.  TN  37388 


Final  Report  for: 

Summer  Research  Program 
Calspan  Corporation/AEDC  Operations 
Arnold  Engineering  Devlopment  Center 
Arnold  Air  Force  Base 
Tullahoma,  TN 


Sponsored  by: 

Air  Force  Office  of  Scientific  Research 
Boiling  Air  Force  Base.  W'ashington,  D.C. 


August  1992 


6-1 


A  CELL  AVERAGED  APPROACH  TO  THE  SOLUTION 
OF  INTEGRAL  CONSERVATION  LAWS 


Blair  H.  Rollin 
Graduate  Research  Assistant 
Department  of  Mathematics 
The  University  of  Tennessee  Space  Institute 


Abstract 

An  analytical  cell  averaging  approach  is  applied  to  the  Local  Lagrangian  Finite  Vol¬ 
ume  method  developed  for  computing  solutions  to  the  compressible  flow  equations.  This 
approach  eliminates  the  need  for  pointwise  evaluation  of  flu.xes  and  coupled  with  nonoscila- 
tory  interpolating  functions  yields  a  highly  accurate,  conservative,  stable  scheme.  This  is 
done  without  the  addition  of  any  terms  not  present  in  the  original  equations,  such  as  arti¬ 
ficial  dissipation  terms.  .Nor  are  the  equations  split  into  characteristic  fields.  The  scheme 
is  derived  and  then  demonstrated  on  two  different  fluid  flow  problems. 


A  CELL  AVERAGED  APPROACH  TO  THE  SOLUTION 


OF  INTEGRAL  CONSERVATION  LAWS 
Blair  H.  Rollin 


Introduction 

The  purpose  of  this  paper  is  to  formally  document  a  new  computational  algorithm  and 
demonstrate  its  viability  for  fluid  flow  applications. 

It  can  be  shown  [1]  that  the  Euler  equations,  describing  the  flow  of  an  inviscid  perfect 
gas,  can  be  written  in  the  form  of  a  homogeneous  coupled  system  of  nonlinear  ordinary 
differential  equations.  From  this  formulation  it  is  seen  that  the  domain  of  dependence 
of  the  system  is  symmetrically  distributed  about  streamlines  of  the  flow  (in  a  differential 
sense).  This  fact  makes  a  Lagrangian  formulation  of  the  equations  a  natural  choice  for  the 
basis  of  numerical  discretization. 

Since  late  1990,  K.  C.  Reddy  of  The  University  of  Tennessee  Space  Institute  and  the 
author  have  been  developing  a  novel  approach  to  solving  the  compressible  flow  equations 
based  on  a  Lagrangian  finite  volume  in  space-time  on  fixed  grids  [2]  [.3].  The  algorithm 
is  particularly  elegant  in  that,  unlike  most  popular  schemes  used  today,  it  includes  no 
additional  terms  not  discretized  from  the  original  equations  to  provide  stability. 

Recently,  by  adopting  a  cell  averaged  approach  to  the  equations,  great  strides  in  accu¬ 
racy  and  stability  have  been  made.  This  paper  focuses  on  the  new  cell  averaged  approach. 


6-3 


I  will  summarize,  in  brief,  the  development  of  the  Local  Lagrangian  Finite  Volume 
(LLFV)  method  since  [2]  was  submitted.  Then  the  present  form  of  the  algorithm  will  be 
derived  and  applied  to  several  test  cases. 


Integrating  V  o  (Q ,  £  )  =  O  over  volume  V  yields  f^y(Q  ,  E  }  o  nds  —  O  ot 


/  -  /  Q(x.t'‘)dx  + 

(.)-l  l) 


F(r, +  1/2(0. n 


".-./2(‘"'^‘)  F(r^_xn{t).t) 


A-./2«'* 


t^O-1/2  =  O 
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n  +  1 


Figure  1:  A  Lagrangian  Finite  Volume  in 


with 


(2) 


where  x  = 


/  \ 
0 


F  =  -Qu  4-  £ 


VPU  / 

T‘j±\ii  are  solutions  to  the  differential  equation 
dx 


n+l  \  _ 


—  =  u.  with  r^±i/2(«  ) 


Backward  integration  in  time  gives 


=  ^3±ll2- 


Since 

dr  -  =  (-7-  +  l)^''^dt  =  (u 

dt 

we  may  rewrite  the  previous  equation  as 


(3) 


/  Q( 


1 


u 

_  n 


space-time. 


^j±l/2- 


2-f  l)^/^dt, 


,t’^)dx  4- 
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/  F(^j+i/2(n,<Wi~  /  F(rj_iy2{n.ndi  =  o. 

At’')  At’') 


it”*') 


From  the  definition  of  r 


ft’'*'  ft"*' 

I  r'dt  =  /  u(r{t),  <)£]!« 

Jt"  Jt" 

rt"*' 

^j±l/2  —  •^J±l/2 

Jt" 


Once  ;his  has  been  integrated. 


d  l‘*]+\/2it) 


pdx  =  0. 


In  this  way  the  Lagrangian  formulation  effectively  transforms  the  system  of  three  PDF’s 
to  that  of  three  ODE’s. 

.approximating  the  integral  on  the  right  by  the  trapezoidal  rule  yields 


(4) 


^j±l/2  =  ^j±l/2  - 


Uj±l/2  —  tlx^j±i/2,t  ))/2. 


With  this  approximation  we  are  assuming  the  characteristics  0±t/2  straight  lines  of 

slope  —  =  Ujj.if2-  Integrating  the  second  two  terms  in  equation  3  similarly,  we  have. 


(5) 


*J+I/2 


f^}  + 


Q{x,  )dx 


12 


„  —  — 

—  Q(x,t  )dx  +  jF" 2+1/2 At  —  Fj_i/2At  =  O  . 

5;-1/20>) 


where 

F2±i/2  =  (F;^;/2+^(^2ii/2,«”))/2. 
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and  the  so  called  flux  vector  F  is  defined  by  equation  2. 


To  preserve  conservative  properties  in  time  and  space  for  discrete  equations  we  choose 
a  solution  space  at  each  time  step  on  a  fixed  spatial  grid  for  approximating  our  solution  to 
the  PDE.  and  use  this  solution  space  consistently  throughout  each  stage  of  the  computation 
for  integral  evaluation  [3]  . 

We  choose  a  solution  space  that  is  piecewise  constant  within  each  interval  ( Xj_i/2,  -c  ,*1/2 ) 
for  integration  purposes  and  use  linear  interpolation  of  quantities  between  nodes.  This 
choice  yields  a  conservative  stable  scheme  which  can  be  shown  to  be  dissipative  of  order  2 
on  the  scalar  equation  +  aux  =  0,a  =  constant.  Also,  this  provides  only  0(h)  accuracy. 

To  increase  accuracy  to  0{h^),  a  piecewise  linear  continuous  solution  space  w?s  chosen 
for  both  integration  and  interpolation  purposes.  This  method  proved  highly  accurate,  but 
unstable  for  some  initial  conditions  on  the  vector  equations.  Following  ideas  of  the  so-called 
MUSCL  and  flux-limited  schemes,  the  previous  two  results  suggested  that  a  piecewise  linear 
slope-limited  approximation  might  provide  stability  along  with  high  accuracy.  However, 
this  method  also  proved  to  be  unstable  for  certain  initial  conditions,  presumably  due  in 
part  to  the  following. 

It  was  observed  that  when  the  piecewise  constant  solution  space  was  used  for  both  in¬ 
tegration  and  interpolation  the  method  was  also  unstable  (recall  that  this  choice  is  stable 
if  linear  interpolation  is  used).  W'hen  a  discontinuous  solution  space  is  used  for  interpo¬ 
lation  an  ambiguity  arises  in  determining  a  function  value  at  a  discontinuity.  Moreover. 
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precise  iocation  of  the  argument  is  critical  in  the  sense  that  a  very  small  change  in  the 
argument  could  produce  a  very  large  change  in  the  corresponding  ordinate.  This  is  clearly 
undesirable  in  an  iterative  process. 

With  these  thoughts  in  mind  we  now  develop  the  present  form  of  the  algorithm. 

Let  X  6  ( -r  -t' .+1/2)  ^ttd  form  a  LLF\'  bounded  by  x  —hj’l  and  x  -  h  Tl  at  .  We 
have 


(6) 


/: 


x+hn 


z—fi/2 


Q(y.t''^')dy  -  J 


iR(i) 


Jrfiit’' 


iK't.Uil 


■(ITfl 


-I 


Q(  I/.  t"')dy  + 

u(^) 

F(rc(f).n 


■dTi  =  O. 


with  Tfl  the  solution  to 


^  =  it  =  z  +  k/2. 

at 


where  we  define 


^r[x)  = 


with  similar  equations  for  Since 


/  'V:  dTn,L^  /  F(ra£,(0, 


we  rewrite  equation  6  and  average  over  cell  j  to  get. 


(7) 


Q{yj'''*')dydx~-  /  Q  { y,  t’')dydx  + 

^3-1/2  ■^r-h/2  h  Jiiiz) 


1  r^t-h/2 

^  ■'^3-1/2 

/  /  F(TR{t].t]dtdx--  /  F{TL{t).t)dtdx  =  O 

dt’'  n  Jt” 
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Using  the  midpoint  rule  we  have  the  following  results  for  the  first  two  terms 


-  /  /  Q(y,t  ^  )dydx  -  /  Q(x.(  )<fi  +  C?(/z‘) 

-  /  /  Q(y,f  )<iydx  =  /  Q(x,t  )di  + 


provided  Q  is  sufficiently  smooth. 


Consider  the  remaining  term 


1  r^]'¥xi2  r*  ^ 

7/  /  F(rfl(<),ndfdx 


and  its  counterpart  with  ti  replacing  tr.  When 

/ 

T-fi(  «"■*■*)  =  2;j-1/2  +  fi/2  =  X; 

X  =  < 

=  X_,_i/2  -  /l/2  =  Xj.i, 

and 

I  =  ^2+1/2  +  hl2  = 

=  3:j+i/2  -  h/2  =  Xj. 

So,  for  example,  the  first  term  signifies  the  double  integral  of  F  over  the  region  in  the 
X  —  t  plane  bounded  in  time  by  the  lines  t  =  <"  and  t  =  and  bounded  spatially  by 
streamlines  intersecting  t  =  at  Xj  and  Xj+i.  Notice  that  this  formulation  eliminates 
the  need  for  pointwise  evaluation  of  fluxes.  As  we  shall  see,  these  approximations  lead  to 
discrete  equations  which  will  be  solved  iteratively. 

Recall  the  result  S]±\I2  -  ^j±ii2  ~^;±i/2^t-  equation  4.  This  result  may  be  achieved 
through  another  interpretation  which  will  prove  useful  in  our  current  pursuit;  namely,  to 
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eliminate  pointwise  evaluation  of  quantities.  We  reinterpret  the  result  in  the  foliowine; 


manner.  Redefining 


1  ! 

^}±i/2  =  ^j±in  ~  r  -  /  “(T;±i/2(r)- ndr-±i/2At. 

sj±i/2  ttiay  be  interpreted  as  Xj±i/2  minus  the  average  value  of  u  along  the  curve  of  inte¬ 
gration  of  the  flux.  fj±i/2  multiplied  by  At.  Then 


>2:1/2  ~  ( 


/(n  I  u*(  O  i"  1)  ^  df 


,n+i 

Jt" 


Ml 


2(n.n(u'  +  n  diAi 


since  c/r  =  (u^  +  1)  ^  dt.  Choosing  some  numerical  approximation  to  u  along  “ 

where  the  subscripts  on  u  are  suppressed  for  convienence. 

.(H+l 

/  (u'(r2±i/2(0.0+  =  !l(d>l)!lAl 

jf' 


ant 


I  u.(Tj±i/2(t).  t)(u'  +  D^^^dtAt  =■  u||(  u,  1 

since  '^  ■=  ii  =■  constant.  The  notation  <5~jii/2  simply  means  the  distance  between  the  two 
points  (^j±i/2)^'‘)  and  ( J2±i/2>  in  tbs  x  —t  plane.  So 

^>±1/2  =  ^]±l/2  —  ~  -^>±1/2  —  “j±l/2‘^^ 

in  agreement  with  equation  4. 

With  this  motivation,  we  define 


^J±l/2  =  -^2*1/2  -  •^]±l/2‘^i- 
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Figure  2:  Region  of  flux  integration 

where  IT  is  taken  to  be  the  average  value  of  u  in  the  region  of  integration  of  the  fluxes 
shown  in  figure  2.  VVe  may  now  formally  define 


/■*J+l/2 

^3  +  1/2  - 

/  dTf{dx. 

f ’•«(<’*) 

r*)-nn 

1 

II 

/  dridx. 

''^7-1/2  ■ 

Jr  tin 

It  is  now  a  simple  extension  to  include  a  source  term  in  our  equation.  For  Q,  £?x  =  S . 
we  have 

d  d  r  .  . 

—Q  +  —^E  -  /  sdr})  =  O, 

Ot  Ox  Jxo 

where  xq  is  some  fixed  point  in  our  domain^  Before  averaging,  our  surface  integral  expan¬ 
sion  will  now  include  the  new  term 

-/  /  S{n,t}dr}dt -i-  /  S(r}J)dT}dt. 

Jt’'  JlQ  Ji" 

'It  may  be  rioted  that  leaving  the  source  term  on  tlie  right  side  of  the  original  equation,  integrating  over 
a  LLFV  centered  at  x.  averaging  this  term  and  applying  the  midpoint  rule  produces  the  same  result. 
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Moving  these  terms  to  the  right  side  of  equation  d,  averaging  the  equation  over  ceil  j  and 
applying  the  midpoint  rule  to  the  resulting  integral  on  the  right  leaves  us  with  equation  <. 


with  a  nontrivial  right  side  given  by 


RHS  =  j^  J  S(x.t)dxdt  - 


(8) 


/2(d  ft’'*' 


S(x,  t)dxdt. 


This  is  simply  the  integral  of  S  over  a  LLFV  centered  at  ;  i.e.,  region  Rj  in  figure  3. 


■^;+i 

"7 - ^ - r~ 

J -I  j  J+l 


n  +  1 


n 


Figure  3:  Region  of  source  integration. 

It  is  now  clear  that  we  will  need  to  integrate  various  dependent  variables  over  the 
regions  shown  in  figures  2  and  3.  For  simplicity  we  consider  the  integration  of  the  generic 
dependent  variable  /  over  the  region  Rj  in  figure  3.  At  any  stage  of  the  iteration  we  have 
an  estimate  of  the  solution  at  time  and  we  know  the  solution  at  Because  of  the 
Lagrangian  formulation  it  is  evident  that  assuming  the  solution  to  be  constant  along  any 
given  r  in  ((n  -  I/2)At.(u  +  l/2)Al)  would  be  a  better  assumption  than  assuming  it 
constant  along  x  —  constant  in  the  same  time  interval.  Therefore,  we  divide  the  rcj,iOrt  /?, 
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into  two  parts  shown  in  figure  3.  Rj,  that  region  of  in  ((n  +  l/2)At,{n  +  1  lAt)  and 
that  region  of  Rj  in  (nAi,(n  4-  l/2)Af).  We  now  make  the  approximations 


^  ^  ^  “  *4-./ 

1  /•’■«(<")  1  fi, 

^  JrJe.. 


/2 

/2 


tj  +  l/j 


fix,t^'*'^)dxRj 

fix.t’^UxR^ 


AT{t)  -  Tfiit)  -  rjr(0 


i.e.;  /  is  assumed  constant  in  the  regions  Rj  and  R^^  with  value  equal  to  the  average  value  of 
/  at  in  ^j+1/2)  1  3,nd  the  average  value  of  /  in  (^^-1/2,  ($'^+1/2)  at  respectively. 

The  integral  of  /  over  the  regions  'Rj±i/2  •  Rj±i/2 's  defined  entirely  similarly. 

The  preceding  discussion  has  now  led  us  to  the  following  approximation  to  equation  7 
including  the  source  terms  given  by  equation  8. 


(9) 


/  -  /  Q(x.r)dx  + 

•^^1-1/2  Mi-xn 


■j+i/ 


-  f’*'  F(x.t'")dx  +  r'""  - 

.'■>  h  Jx,  '  ’ 


i{  ^ 

i;  ^ 

r^j+i/a 


T'  F(x.r)dx  +  r  F(x,r-^^)dx}  ~ 

/ 

^  ^^2  4-. 


ji  h 


/■^7  +  I/3 


5(x,r+Mdx}  =  O 


where  —  ^j+i  —  and  S^j  —  ^j+i/2  —  ^j_i/2-  Note  that  the  R's  can  be  written  in 

terms  of  the  ^’s  as  follows. 


—  ZhAt  At 

Rj+U2  -  ^ 


—  3hAt  At  ^ 

~  +  —(^^+1/2  -  ^7-1/2)- 


6-13 


hilt  3At  hAt  :}Ai 

K}+II2  =  +  l  -  ^j)  E.}  =  -J-  +  -y  (^;  +  l/2  -  ■t,-l/2)- 

Equation  9  is  the  LLFV  method.  We  now  have  everything  sufficiently  defined  so  that  the 
iterat”'e  process  may  begin. 

We  need  only  choose  a  solution  space  on  which  to  perform  the  integration  required  in 
equation  9  and  a  method  of  integration. 

For  the  solution  space  we  choose  interpolating  functions  devised  by  Harten  and  Osherfdj. 
These  are  piecewise  linear  slope-limited  in  each  cell,  constructed  from  the  cell  averages  of 
Q,  and  provide  0(h~)  appro.ximation  to  the  components  of  Q.  Moreover,  they  have  the 
additional  property  that  the  cell  averages  of  the  interpolating  functions  are  an  Oih^)  ap¬ 
proximation  to  the  cell  averages  of  the  components  of  Q.  These  interpolating  functions 
can  be  shown  to  be  nonoscillatory  in  the  sense  that  the  number  of  local  extrema  in  the 
solution  to  the  scalar  equation  will  not  increase. 

Because  of  the  form  of  the  interpolants.  the  first  two  terms  in  equation  9  can  be 
evaluated  exactly. 

The  remaining  terras  are  integrated  by  two  point  Gaussian  quadrature  in  each  subin¬ 
terval  in  which  Q  is  continuous.  This  provides  0{h^)  accuracy  for  sufficiently  smooth  flux 
and  source  functionsfo]. 

As  stated  previously,  because  of  the  coupling  of  the  solution  at  time  level  n  -f  1  and  the 
location  of  the  ^'s  at  time  level  n.  it  is  necessary  to  solve  the  resulting  implicit  equations 
iteratively.  The  Locally  Implicit  Method  (LIM)  developed  by  Reddy  and  Benek[6]  was 
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chosen  as  the  method  of  solution. 

Results 

The  algorithm  was  applied  to  two  test  cases,  the  familiar  shock-lube  problemfl]  and  the 
quasi-one-dimensional  nozzle  flow  problem[7],  both  of  which  have  analytical  solutions. 
Shock  Tube: 

We  wish  to  solve  the  Euler  equations,  equation  1,  on  an  infinite  domain  with  initial 
conditions  of  two  constant  states  separated  by  a  discontinuity  in  the  components  of  Q. 
In  equation  9  F  is  given  by  equation  3  and  S  =  o.  Rather  than  present  the  details.  I 
merely  state  that  sweep  dependent  coefficients  for  the  LIM  were  derived  through  a  lengthy 
analysis.  The  coefficients  used  were 

2i2, 4.1/2  —  1/2 

- Cj  left  to  right  sweeps 

and 

^  ,  I  Rj  +  I/2  -  I 

iix  -I - c,  right  to  left  sweeps 

where  c  is  the  speed  of  sound. 

The  density  solution  for  initial  conditions  pi  =  10^;///;,  =  l;pH  =  10^;///?  =  0.01;  = 

Ufi  =  0,  is  shown  in  figure  4.  Note  that  the  solution  is  highly  accurate  but  virtually 
nonoscillatory.  The  contact  discontinuity  is  particularly  well  resolved. 
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Shock  Tube  Flow;  NGRD=100  CFL=2.0  TIME=3.9msec 


Figure  4;  Density  in  the  shock  tube. 


Nozzle  Flow: 

The  quasi-one-diraensional  nozzle  flow  equations  can  be  written  as  Q(  4-  =  S  where 

(  \ 


Q  =  a 


(  \ 

P 

pu 


E  =  ft 


pu 


pu'  +  p 


.  5  = 


/  \ 
0 


pa 


\  »  / 


y  e  y  (e  +  p)u  y 

and  a{x)  is  the  crossectional  area  of  the  nozzle.  The  equations  were  solved  in  a  diverg¬ 
ing  duct  with  domain  x  €  [0,  lOj  and  the  area  of  the  duct  given  by  a(x)  =  1.398  -f- 
0.347 arctan[0.8(i  —4)].  Empirical  knowlege  has  indicated  that  the  same  coefficients  used 
for  the  shock  tube  problem  are  appropriate  and  work  well  for  the  nozzle  flow  problem. 


These  were  used  and  indeed  did  provide  convergence.  Figure  5  shows  the  steady  state 
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Nozzle  Flow:  Steady  State;  N(j!RD=o() 
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Figure  .5;  Mach  number  in  diverging  duct. 

Mach  number  distribution  for  inlet  conditions  To  =  300A’ ;  po  =  I6ar;  A*  =  0.8  with  outlet 
boundary  conditions  to  locate  a  shock  at  x  =  4.  No  convergence  studies  were  undertaken. 
It  was  necessary  to  use  a  CFL  number  less  than  or  equal  to  1  due  .;o  the  difficulty  in  im¬ 
posing  boundary  conditions  with  large  time  steps,  .\gain  th'  solution  is  highly  accurate. 

Conclusions 

The  purpose  of  this  research  was  to  formalize  the  new  cell  averaged  algorithm  and  inves¬ 
tigate  its  viability  on  some  fluid  flow  applications. 

The  new  cell  averaged  approach  applied  to  the  LLFV  method  provides  a  mathematically 
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elegant,  highly  accurate,  numerically  stable  method  of  solution  for  conservation  laws.  In 
order  for  it  to  be  applied  to  realistic  flow  problems,  future  work  must  concentrate  on 
implementing  boundary  conditions  in  conjunction  with  large  time  steps  and  extensions  to 
multiple  dimensions. 

Acknowledgments 

I  wish  to  thank  .\FOSR  and  RDL  for  flnancia  support  of  this  research.  I  also  wish  to 
express  my  gratitude  to  all  those  at  C.\LSPAN  Corp./.AEDC  Operations  who  have  assisted 
me  over  the  last  two  summers,  in  particular  Stephen  L.  Keeling  and  Robert  \V.  Tramel  for 
insightful  discussions.  Finally  I  wish  to  thank  K.  C.  Reddy  .  without  whose  ideas  none  of 
this  research  would  have  been  possible. 

References 

[1]  Hirsch.  C.,  Numerical  Computation  of  Internal  and  External  Flows.  Volume  3.  John 
Wiley  and  Sons,  West  Sussex,  England,  1990 

[2j  Rollin,  B.  H.,  ,-l  Local  Lagrangian  Model  for  the  Infinite  Domain  Shock  Tube  Problem. 
AFOSR  Summer  Research  Program,  Final  Report.  1991 

[3]  Reddy,  K.  C..  Rollin.  B.  H..  .1  Local  Lagrangian  Method  for  Conservation  Laws. 
Developments  in  Theoretical  and  .\pplied  Mechanics.  Volume  XVI.  Proceedings  of 


6-18 


the  Sixteenth  Southeastern  t.'onference  on  Theoretical  and  Applied  Mechanics.  ed.i. 


B.  N.  Antar,  R.  Engels.  A,  A.  Prinaris,  T.  H.  Moulden.  The  University  of  Tennessee 
Space  Institute.  1992 

[4j  Harten.  .4.,  Osher.  S..  Uniformly  High-Order  Accurate  .Vonoscillatory  Schemes.  SL4.M 
Journal  of  Numerical  Analysis.  V'ol.  24.  No.  2.  April  1987 

[5]  Davis.  P.  J..  Rabinowitz.  P..  Methods  of  Sumerical  Integration.  .\ca.dcmic  Press.  New 
York.  1975 

[6]  Reddy.  K.  C..  Benek.  J.  .4..  .4  Locally  Implicit  Scheme  for  .UD  Compressible  Viscous 
Flows.  .4IAA-90-1525.  June  1990 

[7]  Zucker.  R.  D..  Fundamentals  of  Gas  Dynamics.  Matrix  Publishers.  Chesterland.  OH. 
1977 


6-19 


ANALYSIS  OF  ACOUSTIC  OSCILLATIONS  IN  CAVITIES 
WITH  SPOILER  A'lTACI IMENTS 


Daniel  E.  Schatt 
Master  of  Science  Candidate 
Department  of  Aerospace  Engineering 
University  of  Tennessee  Space  Institute 


Final  Report  for: 

Summer  Research  Program 
Arnold  Engineering  Development  Center 


Sponsored  by: 

Air  Force  Office  of  Scientific  Research 
Bolling  Air  Force  Base,  Washington,  D.C. 


August  1992 


7-1 


ANALYSIS  OF  ACOUSTIC  OSCILLA'nONS  IN  CAVITIES 
WITH  SPOILER  ATTACl  IMEN'I'S 


Daniel  E.  Schatt 
Master  of  Science  Candidate 
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Abstract 

An  attempt  was  made  to  predict  the  amplitude  of  acoustic  oscillations  in  cavities  with  various 
types  of  spoilers  being  used  as  suppression  devices.  A  computer  code  w'as  written  for  this 
purpose.  The  basic  approach  was  to  represent  the  spoiler  as  a  thicker  inititil  boundary  layer,  which 
would  have  the  same  effect  in  suppressing  the  acoustic  oscillations.  With  this  equivalent  boundary 
layer,  each  spoiler  configuration  could  be  assigned  a  certain  drag  coefficient,  which  served  as  the 
primary  input  device  for  the  code.  The  prediction  was  made  over  a  wide  range  of  Mach  numbers, 
from  subsonic  to  supersonic.  The  results  of  the  computer  code  were  compared  with  experimental 
data,  and  also  with  empty  cavity  cases  (no  spoiler). 
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ANALYSIS  OF  ACOUSTIC  OSCILLATIONS  IN  CAVITIFS 
WITH  SPOILER  ArrACHMF.NTS 


Daniel  E.  Schati 


iNTS.o.Ducrm 

Aeroacoustic  oscillations  in  flow  over  cavities  have  been  studied  since  the  1950’s. 

Investigation  of  this  phenomenon  has  great  importance  because  cavities  are  encountered  in  a  wide 
variety  of  applications.  A  cavity  is  defined  as  a  cutout  in  a  surface.  A  typical  example  is  a 
weapons  bay  in  a  bomber  aircraft.  In  my  own  work,  this  application  was  the  primary 
consideration. 

Essentially,  the  oscillations  are  due  to  excitation  of  the  instabilities  of  the  shear  layer  w  hich 
develops  over  the  cavity.  Interaction  of  the  shear  layer  with  the  leading  and  trailing  edges  creates 
resonance  at  certain  frequencies,  which  intensifies  the  response  of  the  shear  layer  at  those 
frequencies.  The  oscillations  occurring  in  the  cavity  (weapons  bay)  can  be  so  intense  that  they 
cause  damage  to  sensitive  instrumentation  in  or  on  the  stores.  In  some  cases,  they  can  even  excite 
the  main  structural  modes  of  the  aircraft.  Therefore,  it  is  very  important  to  be  able  to  predict,  at 
least  approximately,  the  characteristics  of  the  oscillations  (i.e.  frequency  and  amplitude).  In  my 
research  work  in  the  1991  Summer  Research  Program,  I  assisted  in  the  development  of  a  computer 
code  (known  as  the  Cavity  Acoustic  Prediction  Code,  or  CAP  Code)  that  could  give  rough 
predictions  of  the  frequency-amplitude  spectrum  up  to  Mach  1.5. 

The  objective  of  my  current  work  was  to  extend  to  the  code  to  different  flow  configurations. 
Frequently,  suppression  devices  are  used  in  an  attempt  to  reduce  the  intensity  of  the  oscillations. 
These  suppression  devices  are  normally  spoilers  which  are  positioned  upstream,  or  at  the  leading 
edge,  of  the  cavity.  By  interfering  with  the  initial  development  of  the  cavity  ...hear  layer,  the 
magnitude  of  the  oscillations  is  drastically  reduced.  My  objective  was  to  incorporate  into  the 
existing  computer  code  a  method  for  predicting  the  magnitude  of  the  oscillations  for  cases  where 
spoilers  are  used.  Several  different  types  of  spoilers  were  tested  experimentally,  and  these  results 
were  compared  with  the  computer  code  results. 
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Before  discussing  my  work  in  detail,  it  is  necessary  to  brietly  review  the  fundamentals  of  the 
CAP  Code.  The  code  is  based  on  the  principle  that  the  so-called  “edgetone”  frequencies,  first 
introduced  by  Rossiter,  largely  determine  the  characteristics  of  the  oscillations.  These  frequencies 
are  the  frequencies  at  which  vortices  are  shed  from  the  leading  edge  of  the  cavity.  ITesc  vortices 
propagate  downstream  within  the  shear  layer  and  interact  with  the  trailing  edge.  Tliis  generates 
pressure  pulses  which  travel  upstream,  interact  with  the  leading  edge,  and  generate  the  vonices. 
Thus,  we  have  a  continuous,  self-sustaining  feedback  loop.  Rossiter  devised  an  empirical  formula 
based  on  this  -^odel,  to  determine  the  edgetone  frequencies; 


Vo.  (m  -  y) 


where  m  =  1,  2,  3,  ...  =  the  frequency  mode  number  of  the  edgetone. 
tj  iSi  =  empirical  parameters 

a^o  ==  speed  of  sound  based  on  freestream  static  temperature 
aj  =  speed  of  sound  based  on  freestream  total  temperature 


As  can  be  seen,  the  edgetone  frequencies  occur  in  integral  modes,  similar  to  harmonics.  The 
frequencies  observed  in  cavity  oscillations  occur  invariably  on  the  edgctoncs,  with  the  first  three  or 
four  modes  dominating,  and  the  higher  modes  fading  into  the  background  noise. 

The  basic  assumption  of  the  CAP  Code,  justified  by  experimental  observation,  is  that  when 
the  edgetone  frequencies  are  at  or  near  the  natural  acoustic  frequencies  of  the  cavity,  resonance 
occurs  and  the  edgetones  are  greatly  amplified.  Frequencies  away  from  the  edgetones  are  damped 
out,  with  the  degree  of  damping  determined  empirically.  Specifically,  the  damping  is  determined 
by  the  so-called  damping  ratio.  Each  edgetone  frequency,  for  a  given  case,  has  associated  with  it  a 
damping  ratio,  which  is  expressed  as  a  function  of  Mach  number  and  damping  factor.  The 
damping  factor,  m  turn,  depends  on  the  edgetone  frequency.  In  order  to  determine  the  frequency- 
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amplitude  spectrum,  frequencies  through  the  spectrum  tu-e  tested  at  set  interv-als.  I'or  each 
frequency  that  is  tested,  the  amplitude  is  determined  using  the  first  ten  damping  ratio  corresponding 
to  the  first  ten  edgetone  frequencies,  and  the  maximum  is  extracted  as  the  amplitude  for  that 
spectrum  frequency.  The  amplitude  is  determined  by  calculating  a  respon.se  coefficient  R,  and  then 
multiplying  by  a  reference  pressure. 

The  details  of  this  analysis  can  be  found  in  last  year’s  report.  What  1  wish  to  focus  on  now  is 
how  the  initial  boundary  layer  thickness  affects  the  calculated  amplitudes  for  the  various 
frequencies.  The  parameter  0^  in  Rossiter’s  equation  was  expressed  empirically  as: 

04  =  (0.6163  +  0.0178  M»)  (l  -  e  ' 


where  >jp  =  o  — ,  a  turbulent  mixing  position  parameter 


and  a  =  the  similarity  parameter  for  turbulent  mixing,  after  Bauer 

Therefore,  the  boundary  layer  thickness  affects  the  edgetone  frequencies,  which  in  turn  affects 
the  amplitude  (sound  pressure  level). 

This  fact  was  used  in  the  analysts  of  the  spoiler  cases.  First,  the  initial  boundary  layer 
thickness  approaching  the  cavity  is  assumed  to  be  the  turbulent  boundary  layer  thickness  according 
to  the  results  of  Whitfield  and  Tucker.  This  thickness  is  a  function  of  Mach  number.  These  values 
are  inputted  into  the  CAP  Code,  which  computes  the  sound  pressure  levels  of  the  empty  cavity 
case  (no  spoiler)  for  various  Mach  numbers.  It  is  important  to  point  out  that  here,  and  henceforth, 
we  arc  referring  to  the  overall  sound  pressure  level,  which  is  essentially  a  root-mean-square  of  all 
the  amplitudes  corresponding  to  all  the  frequencies,  for  a  particular  case.  These  computed  sound 
pressure  levels  (SPLs)  are  compared  with  the  experimental  values,  for  each  Mach  number.  This 
allowed  the  difference  between  the  CAP  Code  and  experimental  values  to  be  calculated,  for  each 
Mach  number. 

These  empty  cavity  results  were  then  used  for  the  spoiler  cases.  The  differences  in  sound 
pressure  level,  for  the  various  Mach  numbers,  were  subtracted  from  the  experimental  spoiler 
results  to  obtain  an  equivalent  SPL  for  the  CAP  Code.  Then,  boundary  layer  thicknesses  were 
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inputted  in  an  iterative  fashion  until  the  CAP  Code  SPL  matched  this  equivalent  SPL.  'I'his 
boundary  layer  thickness  was  taken  as  an  equivalent  thickness;  in  other  words,  the  spoiler  can  be 
considered  to  have  the  same  effect  on  the  SPL  as  a  greater  initial  boundtiry  layer  thickness,  as 
calculated  by  the  CAP  Code. 

It  is  appropriate  at  this  stage  to  briefly  review  the  various  spoiler  types  used  in  the  analysis 
(and  compared  with  experiment).  Two  spoilers  with  sawtooth  edges  were  used:  one  with  coarser 
sawtooth  and  one  with  finer  sawtooth.  A  solid  spoiler  and  porous  spoiler  were  also  used.  Finally, 
a  flap  configuration  was  tried  in  two  different  sireamwise  positions.  All  of  these  configurations  arc 
illustrated  below; 


CONFIGURATION  120 
ja  SAIMTOOTH  SrOM.II(. 
Xll  mO 


CONFIGURATION  112 
JS  SOHO  SFONlI*. 

XII  m  0 


CONFIGURATION  1M 
M  FOIIOUS  SFOAIR.  (M  %) 
X/l  •  0 


CONFIGURATION  122 

}»  fLAF  SFON.1II. 

Kll  m  0 
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It  will  be  seen  shortly  that  height  of  the  spoiler  is  an  important  parameter.  All  configurations  have 
the  same  height,  although  the  sawtooth  spoilers  and  the  porous  spoiler  used  an  equivalent  height 
due  to  the  geometry. 

An  important  objective  was  to  obtain  equivalent  drag  coefficients  for  the  various  spoiler  types, 
which  is  the  drag  coefficient  resulting  from  replacing  the  spoiler  with  an  equivalent  initial  boundary 
layer  thickness  which  produces  the  same  sound  pressure  level.  In  order  to  do  this,  an  equation  for 
drag  coefficient  was  derived  by  Calspan  engineer  Robert  Bauer,  which  expressed  the  drag 
coefficient  in  terms  of  upstream  and  downstream  boundary  layer  thicknesses;  then  the  downstream 
thickness  is  considered  to  be  the  equivalent  initial  thickness.  This  derivation  was  performed  in  the 
following  manner: 
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These  drag  coefficients  are  based  on  the  dynamic  pressure  of  the  freestream,  q^,  according  to: 

Since  q  =  0.5  f  U^,  the  q  is  obviously  changing  over  the  height  of  the  spoiler.  In  order  to  calculate 
a  more  accurate  Cu,  it  is  necessary  to  correlate  it  to  an  average  q  ( q )  over  the  height  of  the  spoiler. 
Once  q  is  known,  the  new  adjusted  €□  can  be  determined  because  drag  remains  constant  and, 
therefore,  the  product  Cq  q  remains  constant.  Thus: 
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The  difficulty  arises  in  calculating  an  average  q.  This  was  accomplished  by  integrating 
q  =  0.5PU2  over  the  height  of  the  spoiler  using  a  typical  Xfl  power  turbulent  boudary  layer  profile, 
and  dividing  by  the  spoiler  height.  This  procedure  resulted  in  an  analytical  expression  for  the  ratio 
of  q  to  q^.  In  turn,  these  results  were  used  to  adjust  the  drag  coefficients  for  all  the  cases 
computed  previously.  The  derivation,  carried  out  with  the  help  of  Blair  Rollin,  another  .Summer 
Program  participant,  is  shown  below; 
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To  my  knowledge,  this  is  the  first  time  an  analytical  solution  has  been  presented  for  the  total 
momentum  in  a  turbulent  boundary  layer.  In  order  to  verify  this  analytic  solution,  the  results  were 
compared  to  the  results  obtained  from  the  following  derivation  by  Bob  Bauer: 


k  -  S  •  '■j  -  j  ^  ^  ^ 

1 


2  U 


TTie  parameters  ^  ^  can  be  found  in  standard  tables  for  various  1/n  power  profiles. 
An  example  comparing  results  of  the  two  methods  is  shown  below: 
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This  approach  was  integrated  into  the  CAP  Code.  For  each  case,  a  C,,  is  inputted  using  the 
previous  results  as  a  guide  (Co  based  on  q^.  In  addition,  the  turbulent  boundary  layer  thickness 
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;md  the  height  of  the  spoiler  arc  inputted.  The  sawtooth  .spoilers  and  the  porous  spoilers  are  given 
equivalent  heights  (total  area  divided  by  width).  Then,  an  equivalent  downstream  boundary  layer 
thickness  is  calculated.  This  formula  was  derived  in  the  following  manner; 
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As  can  be  seen,  the  is  adjusted  according  to  the  aforementioned  correction.  Using  this 
boundary  layer  thickness,  the  CAP  Code  goes  on  to  calculate  the  overall  SPL. 

These  results  arc,  of  course,  highly  dependent  on  the  inputted  Cq.  It  was  decided  that  the 
difference  in  sound  pressure  level  between  the  empty  cavity  case  and  the  spoiler  case,  rather  than 
the  absolute  sound  pressure  level,  was  of  primary  concern.  Therefore,  the  change  in  sound 
pressure  level  was  plotted  as  a  function  of  Mach  number.  This  was  done  for  both  the  experimental 
and  CAP  Code  cases.  For  each  spoiler  type,  the  CAP  Code  results  were  calculated  for  several 
different  drag  coefficient  values.  The  plots  are  shown  on  the  following  pages. 


LVL 


FLAP 


The  following  conclusions  cun  be  drawn: 

1)  From  the  preceding  plots,  it  can  be  seen  that  the  frecsiream  C„  v\  hich  is  inputted  should  be 
higher  than  0.2  for  lower  Mach  numbers  and  low'er  than  0.2  for  higher  Mach  numbers. 

2)  The  suppression  effectiveness  is  clearly  a  strong  function  of  Mach  number,  and  therefore 
the  equivalent  drag  coefficient  must  change  with  Mach  number. 

3)  With  a  sufficient  experimental  data  base,  the  correct  Cj^s  to  input  could  be  estimated  more 
precisely. 
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Numerical  Modelling  of  Mixing  and  Reacting  Flowficlds 


Paul  Vitt 
Graduate  Student 

Department  of  Mechanical  and  Aerospace  Engineering  and  Engineering  Mechanics 

Universitv'  of  Missouri  -  RoUa 

Abstract 


The  computational  fluid  dynamics  (CFD)  modelling  of  complex  mixing  and  reacting  flow-fields  is  the  goal 
of  the  current  research.  The  mixing  and  reacting  modelling  has  been  broken  up  into  se^-cral  parts  to 
evaluate  their  individual  contributions  to  the  solution.  This  fits  into  the  general  framework  of  ex  aluating 
the  computational  code  GASP  for  engineering  design  purposes.  The  first  part  of  the  study  involves 
qualitatively  evaluating  the  effect  of  turbulence  on  chemical  kinetics  within  the  numerical  modelling.  A 
brief  look  at  the  effects  of  turbulent  Schmidt  number  is  also  presented.  The  selection  of  an  appropriate 
chemistry  kinetics  model  is  very  important  for  flowfletds  where  capturing  the  flame  and  ignition  point  are 
important,  as  is  shown  through  a  comparison  of  shock -induced  combustion  numerical  e.xperiments.  The 
code  is  also  used  to  predict  the  flowfield  of  premixed  hydrogen  air  burner,  which  involves  subsonic  flow 
over  a  back  step  with  an  ignition  torch.  This  case  is  to  provide  a  comparison  with  another  code  which  is 
modelling  the  same  problem.  The  mixing  part  of  the  physics  modelling  is  addressed  through  two  low 
angle  wall  jet  injectors,  in  which  the  GASP  predictions  were  compared  with  c-xperimental  and  other  CFD 
results.  The  main  conclusion  from  this  part  of  the  research  is  that  the  accuracy  of  the  numerical 
simulation  in  GASP  needs  to  be  improved  before  more  complex  flowflelds  can  be  modelled  with 
confidence.  The  chemistry  turbulence  interaction  needs  further  attention,  as  does  the  selection  of  an 
efficient,  accurate  chemistry  model.  Finally,  the  turbulent  diffusion  model  needs  to  be  tuned  before  the 
modelling  of  swept  ramp  scamjet  injectors  (the  final  part  of  this  research)  is  undertaken. 


8-2 


1.  INTRODiyCTTON 

The  development  of  computational  fluid  dynamics  (CFD)  as  a  reliable  tool  for  engineering 
analysis  of  aerodynamic  and  propulsive  flowficlds  requires  that  the  codes  be  validated  against  known 
physical  results,  i.e.  experimental  data.  While  the  aerodynamic  part  of  CFD  modelling  is  fairly  well 
represented,  the  modelling  of  propulsive  flowficlds  (mixing  and  reacting)  is  still  an  area  of  uncertainty. 
To  address  this  problem  with  reprd  to  a  single  code,  the  (jeneral  Aerodynamic  Simulation  Program,  or 
GASP,  several  test  cases  have  been  dev'cloped.  The  final  goal  of  the  research  is  to  determine  if  the  code 
accurately  models  reacting  flowficlds  well  enough  to  be  used  as  an  analysis  tool  for  designing  the  National 
Aerospace  Plane  (NASP)  combustor.  The  design  of  the  combustor  is  very  important:  with  flight  Mach 
numbers  of  8  to  20,  a  supersonic  combustor  ramjet  (scramjet)  engine  is  a  likely  candidate.  This  engine  is 
airframe  integrated  so  that  the  fordxtdy  of  the  aircraft  provides  compression,  as  well  as  the  cowl  and 
injeaor  struts.  In  order  to  minimize  total  pressure  losses  the  combustion  occurs  supersonically  (this  also 
avoids  the  high  static  temperatures  that  would  be  associated  \vith  reducing  the  flow  to  subsonic  conditions 
for  combustion).  The  flow  velocity  varies  only  a  few  percent  through  the  engine  [1].  but  it  is  the  high 
mass  flow  rate  that  allows  only  slight  changes  in  velocity  to  generate  thrust.  Another  consideration  is  that 
the  combustor  length  must  be  minimized:  1.  to  minimize  the  extreme  frictional  losses  at  high  flight  Mach 
number  velocities,  and  2.  to  reduce  structural  weight  and  cooling  requirements.  All  of  these  conditions 
lead  to  the  faa  that  there  is  only  a  very  short  flow  residence  time  in  the  combustor  (on  the  order  of 
milliseconds):  hence  mixing  and  combusting  the  fuel*  quickly  is  a  necessity.  Since  current  ground 
facilities  are  limited  in  their  range  of  hypersonic  applicability,  CFD  modelling  is  important  in  the  design 
of  the  combustor  for  the  NASP. 

This  brings  the  discussion  back  to  the  reliability  and  accuracy  of  computational  models.  The  test 
cases  that  have  been  analyzed  in  this  research  involve  both  mixing  and  reacting,  for  numerical 
experimentation  as  well  as  comparison  with  data.  The  first  case  examines  the  effect  of  turbulence 
modelling  on  reaction.  A  10  degree  compression  ramp/expansion  is  used  for  this  numerical  experiment. 
The  second  test  case  looks  at  the  Burrows  and  Kurkov  test  case,  and  the  effect  of  turbulent  Schmidt^ 
number  on  the  modelling  of  that  problem.  The  third  case  looks  at  the  computations  of  shock-induced 
combustion  over  a  10  degree  compression  ramp  for  different  chemistry  models.  A  fourth  test  case  looks  at 
subsonic  combustion  of  premixed  hydrogen.  The  fifth  case  looks  at  supersonic  low  angle  helium  injection 
into  hypersonic  flow,  with  the  goal  of  validating  the  mixing  predicted  by  GASP. 

7. 7  Mathematical  Background 

GASP  is  an  implicit  solver  of  the  finite-volume  form  of  the  Navier-Stokes  equations. 
Approximate  factorization  routines  solve  either  full  Navier-Stokes  (FNS)  or  the  parabolized  Navier-Stokes 
(PNS)  form  for  space  marching  through  largely  supersonic  domains.  The  chemistry  is  based  on  empirical 

^probiMv  hydrogen,  becaiue  of  iw  high  specific  ittipusle  (tow  motecuiar  weight) 

^Sc(tinb)  conlrob  tuihuleni  nun  difltuian 
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Arrhenius  rates,  with  the  properties  supplied  form  the  JANNAF  (Gordon-McBridc)  curv’e  fits.  Turbulence 
is  supplied  from  the  Baldwin-Lomax  algd>raic  model.  f3].  or  from  any  of  three  two-equation  models.  The 
chemistry  models  chosen  for  study  here  are  outlined  in  the  Appendix  in  terms  of  species  and  reactions. 
Due  to  space  limitations,  for  further  numerical  discussion  of  GASP,  the  reader  is  referred  to  Ref.  [2]  for 
the  complete  system  details. 

2.  DISCUSSION  OF  RESULTS 

2. 1  Effect  of  Turbulence  on  Chemical  Reactions 

Turbulence-chemistry  interactions  ate  very  important,  since  turbulence  can  control  chemistry 
through  mass  diffusion,  and  chemistry  should  influence  turbulence  through  flowfield  gradients.  In  order 
to  model  complex  problems,  it  i$  necessary  to  make  the  calculations  with  as  few  variables  (equations)  as 
possible,  in  order  to  minimize  computational  requirements.  This  numerical  experiment  is  to  determine 
whether  the  simpler  Baldwin-Lomax  algebraic  turbulence  model  can  interact  with  a  chemistry  model  in  a 
similar  fashion  to  the  2-equation  k-e  model,  and  try  to  qualify  the  effect  of  turbulence  modelling  on 
chemical  kinetics.  The  problem  is  a  10°  ramp  compression-expansion,  with  the  inflow  premixed 
stoichiometric  hydrogen-air.  The  conditions  at  the  inflow  are  shown  in  table  1 . 


Mach  No. 

T(K) 

^wall 

HHRIfflIHI 

6.0 

1573.0 

2000.0 

0.0386 

Table  1:  Inflow  conditions  for  reacting  ramp-expansion  problem 


These  conditions  provided  for  shock-induced  reaction.  The  results  are  summarized  in  the  sketch 
presented  in  Figure  1.  The  locations  of  the  maximum  gradient  in  H2O  mass  fractions  are  drawn  on  the 

flowfield  geometry.  Alt  three  cases  used  Dnunmond  chemistry,  as  described  in  the  Appendix.  The  first 
case  used  Baldwin-Lomax  algd)raic  turbulence  (which  is  very  efficient  computationally),  the  second  and 
third  used  Chien's  model  low  Reynolds  number  k-e  turbulence  (both  cases  add  two  new  partial  differential 

transport  equations  to  the  set  which  must  be  solved.  Ref  (4]).  The  difference  between  the  second  and 
third  cases  is  the  amount  of  free-stream  turbulence:  the  second  case  had  O.OZUjuflQ^,  and  case  3  had 

0.(X)2Uj„flQ^.  The  results  are  reflected  in  the  flame  speed  (the  location  where  reaction  begins,  which  was 

taken  here  to  be  the  maximum  gradient  in  water  mass  fraction).  The  flame  speed  is  driven  by  the  mass 
diffusion  rate,  which  is  controlled  by  the  turbulent  viscosity  through  the  Schmidt  number.  The  flame 

would  initiate  at  the  high  temperature  wall  and  propagate  out  into  the  flow  at  the  mass  diffiision  rate  (of 
OH,  the  progenitor  of  H2O)  until  the  temperature  falls  and  reduces  the  reaction  rates,  halting  the 

production  of  OH.  The  algebraic  turbulence  underpredicted  the  flame  speed,  as  can  be  seen  by  the  flame 
front  being  closer  to  the  lower  wall.  In  this  case  the  flame  was  limited  by  the  mass  diffusion  rate  of  OH. 
The  high  freestream  turbulence  case  (2.0%  case)  had  a  flame  speed  that  was  faster  than  the  flow  velocity, 
and  was  limited  only  by  the  ignition  source  which  was  the  shock  wave.  The  third  (low  intensity;  0.2%) 
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case  showed  a  flame  speed  that  was  located  about  half  way  between  the  first  two  cases,  again  being 
limited  by  the  mass  diffusion  rate  of  OH.  It  can  be  seen  that  the  freestream  turbulence  has  a  large 
influence  on  chemistry,  and  is  a  modelling  problem  for  CFD  simulations  currently.  The  question  is:  how 

much  freestream  turbulence  is  there?  Possibly,  to  improve  the  validit>’  of  the  algdtraic  solution,  if  the 
amount  of  freestream  turbulence  is  known,  a  constant  should  be  imposed  on  the  flowflcld.  simulating 

the  freestream  turbulence  of  the  k-e  model. 

2.2  Influence  of  Turbulent  Schmidt  Number 

In  order  to  qualify  the  effect  of  turbulent  Schmidt  number,  which  relates  the  turbulent 
momentum  difiusion  rate  to  the  mass  diffusion  rate,  the  Burrows  and  Kurkov  supersonic  tangential  H2 

injection  case  was  modelled  (Figure  2a.b,  [5]).  This  case  was  chosen  because  there  was  good 
experimental  data  on  the  chemistry  compositions  at  the  final  plane,  and  it  had  been  previously  modelled 
with  GASP,  [2].  The  case  was  first  modelled  with  a  turbulent  Schmidt  number  of  0.7.  which  was  the 
same  as  the  published  results,  and  came  up  with  an  identical  solution  (Figure  2c.  square  point).  Sc^  was 

then  changed  to  0.5,  and  the  case  was  rerun.  The  increased  mass  diffusion  rate  moved  the  flame  front  out 
slightly,  but  not  enough  to  match  the  data.  Tlie  reason  for  the  flame  front  being  closer  to  the  wall  than  the 
experiment  is  that  the  incoming  thick  boundary  layer  is  not  modelled,  and  modelling  this  should  increase 
the  penetration  of  the  mixing  layer  and  flame  front.  The  numerical  experiment  did  qualify  the  amount  of 
effect  that  the  Sc^  has  on  the  kinetics. 

2.3  Influence  of  Modelling  on  Chemical  Kinetics 

The  iKXt  stage  in  the  investigation  is  the  effect  of  the  chemistry  model  on  the  results.  Again,  in 
terms  of  reducing  the  computational  requirement  it  is  desired  to  use  the  least  number  of  chemical  species 
to  capture  the  physics.  The  four  chemistry  models  (detailed  in  the  Appendix)  are  denoted  by; 


Chemistry  Model 

Number  of  Chemical  Species 

Number  of  Chemical  Reactions 

1.  Drummond 

7 

7 

2.  Evans  &  Schexnayder  1 

7 

8 

3.  Evans  &  Schexnayder  2 

12 

25 

4.  NASP  4  (H2/NOjj  extension) 

12 

24 

Table  2:  Denotation  of  chemistry  models 

The  problem  is  a  10®  degree  compression  ramp^,  with  the  inflow  conditions  shown  in  table  3. 


Mach  No. 

^wall 

mam 

HSSBES 

U(m/s) 

6.0 

1273.0 

2000.0 

19600.0 

0,0386 

4934.7 

Table  3:  Inflow  conditions  for  the  reacting  ramp  problem 


^similar  to  the  fint  put  of  section  2. 1 
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These  conditions  were  selected  to  produce  shock-ignition  of  the  premixed  stoichiometric  H2-Air  mixture. 

The  chemistry  models  had  a  larger  variation  in  the  results  than  was  expected  or  desired  (Figures  3.  4.  5). 
The  shock  locations  (as  representative  of  the  flowficld  kinematics)  are  shown  in  Figure  3.  a.-d.  The 
Drummond  chemistry  model  caused  the  shock  location  to  be  located  further  up  the  ramp  but  the  other 
three  cases  were  in  good  agreement  for  shock  location.  Figure  4  shows  water  mass  fraction  contours 
(representing  the  finite  rate  chemistry  kinematics)  for  the  different  chemistry  models.  This  is  where  the 
greatest  disparity  is  shown  between  the  models;  the  flame  speed  must  be  strongly  influenced  by  the 
inclusion  of  the  intermediary  species  in  the  larger  models.  This  case  is  very  close  to  the  ignition  point  of 
the  stoichiometric  mixture,  and  the  smaller  models  (7  species)  are  overpredicting  the  ignition  delay  in  the 
flowfield.  Line  plots  of  the  flowfield  kinematics  are  shown  in  Figure  5  a.  (velocity  profiles)  and  b. 
(pressure  profiles)  at  the  end  of  the  domain,  which  was  intended  to  capture  a  cross-section  of  the  flame 
front  in  the  shock  wave.  OH  mass  fraction  profiles  (Figure  5c)  and  water  mass  fraction  profiles  (Figure 
5d)  at  the  same  location  portray  the  chemical  kinetics.  The  models  all  have  similar  velocity  profiles 
across  the  reacting  shock  layer  (Figure  5a),  indicating  that  the  momentum  physics  is  at  least  common  to 
all  of  the  models.  The  differing  shock  locations  can  be  seen  clearly  in  the  pressure  spikes  (Figure  5b), 
with  the  last  three  models  all  showing  a  similar  location.  It  is  of  note  that  the  first  Evans  &  Schexnayder 
model  (7  species)  predicts  a  shock  location  and  strength  that  is  halfway  between  the  two  larger  (12 

species)  models.  In  terms  of  the  chemical  kinetics,  the  Evans  &  Schexnayder  models  predicted  more 
reaction  (more  OH  and  H2O  production).  The  Drummond  model  and  the  NASP  4  model  were  in  good 

agreement  in  terms  of  the  amount  of  reaction  that  was  present.  The  major  difference  between  the  models 
is  the  ignition  delay,  which  is  evidenced  by  the  larger  models  by  the  reaction  (OH  and  H2O)  extending  all 

of  the  way  across  the  domain,  whereas  the  7  species  models  were  limited  by  the  shock  location.  One 
possibility  for  this  discrepancy  is  the  reactions  rates,  but  the  Evans  &  Schexnayder  models  use  the  same 

reaction  rates  for  the  basic  8  reactions  in  the  7  species  model.  Another  interesting  point  is  that  the  OH 
contours  (not  shown  here)  do  not  show  a  dramatic  increase  in  at  the  start  of  reaction  -  note  in 

Figure  4  c.-d.  that  the  increasing  water  contours  in  the  freestream  are  spread  out  axially^.  This  is 

indicative  of  there  being  only  a  weak  flame  front  in  this  case.  The  presence  of  some  of  the  intermediary 
species  opening  up  new  reactionary  pathways  (with  perhaps  lower  activation  energies)  for  OH  and  H2O 

production  in  the  12  species  models  might  be  the  cause  of  the  increased  reaction  before  the  shock.  The 
conclusion  that  can  be  drawn  from  this  experiment  is  that  the  lower  order  models  do  not  pick  up  the 
ignition  delay  well  in  comparison  to  12  species  model  results,  but  from  the  compositions  after  the 
reaction,  the  7  species  models  do  calculate  a  correct  equilibrium  point. 


^tnoOter  note  of  interest  for  the  kinelia  is  that  the  water  maai  faction  goea  down  inside  the  shock,  as  some  of  the  water  is  turned  back  into 
OH,  and  then  letums  to  its  freestream  equilibrium  value  after  the  shock. 
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2.4  The  Hydrogen-Torch  Problem 

Another  test  case  for  the  chemical  reactions  in  GASP  was  designed  around  a  hydrogen  torch. 
The  torch  was  an  ignition  source  placed  at  the  base  of  an  axisymmetric  backstep.  and  was  to  be  used  to 
light  a  premixed  fuel-air  mix.  The  problem  was  broken  into  three  parts;  1 .  the  air-only  flowTield  was 
established.  2.  the  hydrogen-air  torch  at  the  base  of  the  backstep  was  established,  and  3.  the  inflow  was 
changed  to  a  stoichiometric  hydrogen-air  mixture.  Results  for  the  second  stage  are  shown  in  Figure  6a.-b. 
Temperature  contours  (a)  show  the  limited  influence  of  the  torch  (the  torch  exit  was  about  1/4  of  the 
backstep,  and  the  velocity  vectors  plot  (b)  shows  why;  the  torch  is  buried  at  the  bottom  of  the  recirculation 
zone  after  the  step,  and  the  back  flow  is  enogh  to  confine  the  torch  to  the  area  adjacent  to  the  stq).  Figure 
7a.-c.  shows  the  results  after  ignition  of  the  premixed  hydrogen-air  inflow.  The  velocity  vector/streamline 
plot  (7a.)  shows  that  most  of  the  infow  is  directed  out  through  the  unconfined  top  boundary.  This  is  the 
expected  result;  earlier  cases  with  the  top  boundary  as  a  fixed  wall  thermally  choked  the  flow  at  the 
inflow.  Figure  7b.  and  7c.  show  water  mass  fractions  and  temperature  contours  to  be  increasing  at  about 
the  same  location.  The  ignition  has  spread  from  the  torch  upwards  towards  the  unconfined  boundary. 
From  the  streamlines,  only  about  25%  of  the  inflow  actually  proceeds  into  the  combustion  region  shown  in 
the  plots;  the  rest  is  pushed  out  of  the  top  of  the  domain  by  the  pressure  rise  due  to  reaction.  This  result 
agrees  with  the  expected  physics  for  this  problem,  and  is  demonstrative  of  the  ability  of  GASP  to  model 
these  problems. 

2.5  Supersonic  Helium  Injection  into  a  Supersonic  Stream  (Mixing  Problem) 

Having  addressed,  qualitatively,  the  reacting  flowfield,  the  other  part  of  the  computational  model 
is  correctly  predicting  the  mixing.  Mixing  is  driven  by  the  modelling  of  the  turbulence  in  two  ways;  1. 
directly  through  mass  diffusion  and  the  Schmidt  number,  and  2.  indirectly  through  the  dissipation  of  any 
large  structures  in  the  flowfield,  which  induce  large  scale  mixing  through  warping  of  the  interface.  In 
order  to  evaluate  GASP,  and  its  ability  to  predict  flowfield  mixing,  a  case  with  supersonic  helium 
injection  into  a  Mach  6  flowfield  was  chosen.  This  problem  was  selected  due  to  the  available 
experimental  data,  and  the  flowfield  should  approximate  the  conditions  in  a  scramjet  combustor.  The 
helium  was  injected  through  a  flush  wall  port  at  an  angle  of  IS*’  to  the  plate.  The  computational  mottel 
was  initially  broken  up  into  three  domains;  1.  a  flat  plate  entry  length,  2.  a  nearfield  injection  region,  and 
3.  the  farfield  mixing  region.  Experimental  data  is  from  Ref  [6],  and  GASP  has  previously  modelled  the 
problem  in  Ref  [7].  The  earlier  GASP  results  were  not  as  accurate  as  they  should  be  to  confidently 
undertake  more  complex  flowfield  modelling,  so  this  test  program  was  undertaken  to  see  if  the  accuracy 
could  be  improved. 

Two  separate  cases  were  considered;  an  overpressurized  injection  case  and  a  matched  pressure 
case.  These  characterizations  were  based  on  the  effective  back  pressure  method  of  Schetz  and  Billig  in 
Ref  [8],  as  modified  by  Fuller  to  angled  injection  in  |6].  The  effective  back  pressure  is  simply  an  estimate 
of  the  pressure  that  the  jet  will  see  in  the  flowfield;  Fuller  suggested  that  for  this  case  a  15**  cone  be  used 
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to  rq)resent  the  injection,  and  so  the  effective  back  pressure  is  based  on  the  pressure  after  the  shock  over 
this  cone.  This  is  simply  calculated  using  the  cone-shock  chans.  The  first  case  is  that  of  matched 
injection  pressure,  and  the  second  is  for  5  times  overpressure.  The  inflow  and  jet  conditions  are  shown  in 
table  5. 


T  (K) 

V  (m/s) 

Freestream 

4364.0 

63.3 

957.0 

Matched  P  Jet 

21380.0 

150.0 

1225.5  ^  15® 

HEBSH 

5XPJet 

106898.0 

150.0 

1225.5  ®  15® 

0.3430  (He) 

Table  S:  Jet  and  freestream  conditions  for  the  IS  degree  helium  wail  jet 
The  turbulent  Schmidt  number  was  initially  set  to  0.5,  and  the  algebraic  turbulence  model  of  Baldwin- 
Lomax  was  used.  The  Baldwin-Lomax  model  will  not  accurately  capture  nearfield  turbulent  structures, 
but  by  using  a  wake  modification  through  the  jet  it  is  hoped  that  the  model  will  pick  up  the  downstream 
turbulent  mixing. 

The  results  for  this  initial  investigation  with  GASP  proved  to  not  be  very  accurate.  Figure  8a.-b. 
is  an  axial  plot  of  helium  mass  fraction,  and  shows  the  differences  between  the  IX  matched  pressure  case 
(a)  and  the  5X  overpressure  case  (b).  The  core  penetration  can  be  seen  to  be  much  more  for  the  5X  case, 
as  would  be  expeaed.  In  both  the  experiment  and  the  computations,  the  core  of  the  IX  jet  stayed  down 
along  the  wall  for  most  of  the  domain.  The  mixing  results  are  shown  for  aossflow  planes  of  helium  mass 
fraction.  There  were  experimental  measurements  taken  at  four  stations;  X/D  =  20,  40,  60,  and  80.  The 
first  three  planes  for  the  overpressure  case  are  shown  in  Figure  9a. -c.,  for  the  experiment,  the  curreut 
computation  (labelled  GASP),  and  the  previous  computation  with  GASP,  labelled  FULLER,  et.aP, 
respectively  left  to  right.  At  X/D  =  20.  the  experimental  data  is  very  asymmetric  (otUy  the  left  side  of  the 
data  is  shown),  but  seems  to  indicate  that  the  core  of  the  jet  has  been  split  by  the  overspill  vortices. 
Neither  of  the  CFD  solutions  display  this  effect  and  show  considerably  more  diftusion  and  penetration.  At 
X/D  =  40  and  X/D  =  60,  the  experimental  core  has  rejoined  into  a  single  core,  and  the  CFD  solutions  are 
still  over-difitising  the  jet.  The  Fuller  solution  is  more  diffusive  than  the  current  GASP  solution.  Figure 
10  shows  the  helium  crossflow  contours  for  the  matched  pressure  case;  experiment  and  current 
computations.  The  CFD  core  is  a  little  less  transversely  stretched,  but  overall  the  agreement  is  fairly 
good.  Figure  1 1  shows  the  X/D  =  80  plane  comparing  the  current  CFD  solution  to  the  experimental 
contours.  The  core  is  well  over-penetrated  and  over-diffused.  For  other  types  of  CFD,  a  similar 
experiment  to  CFD  comparison  is  shown  for  the  finite-difference  code  SPARK  [91.  and  also  for  the  Fuller 
solution  in  Figure  12a.-b.,  respeaively.  The  SPARK  solution  has  fair  agreement  with  the  data,  but  the 
Fuller  solution  is  as  far  off  as  the  current  solution.  Figure  13  quantifies  the  observations  above.  Figure 

^doomenled  in  Ref.  |7]  ts  the  refined  grid  aolution 
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Ha  -b  are  conscnauoti  of  mass  plots  (used  for  convergence  enteria).  The  current  domain  is  unconfincd 
in  the  vertical  direction:  it  is  the  bou  shock  wave  turning  the  flow  out  of  the  domain  that  causes  the  large 
axial  drop  in  total  mass  flux  for  both  eases.  The  helium  mass  flux  agreement  is  good  for  both 
overpressure  and  matched  pressure  cases  Figure  He.  shows  the  axial  decay  of  the  maximum 
concentrauon  of  helium  mass  fraction,  comparing  results  for  GASP  CFD  and  experiment.  The  CFD 
solution  ovcrprcdicts  the  decay  rate,  but  this  has  been  documented  for  this  case  for  SPARK  in  Ref,  [9] 
Note  that  both  the  overpressure  and  matched  pressure  CFD  solutions  decay  at  the  same  rate,  whereas  the 
experimental  decay  is  slightly  faster  for  the  matched  pressure  case.  This  is  probably  because  the 
kinemauc  farfietd  is  very  similar  for  the  two  cases  computationally,  and  since  the  turbulent  diffusion  is 
based  on  mean  kinematic  properties,  the  diffusion  for  both  cases  should  be  similar.  This  kinematic 

representation  is  probably  the  cause  of  some  of  the  discrepancy  between  the  solutions  and  the 
e.xpenmen(al  data.  Figure  1 3d.  shows  the  penetration  of  with  a.xial  distance.  The  matched 

pressure  case  is  in  very  good  agreement  with  the  data,  in  terms  of  the  flowTicld  accurately  predicting  the 
core  remaining  along  the  wall  until  some  dow-nstream  location  before  being  lifted  off.  The  penetration  of 
the  overpressure  solution  is  grossly  overpredicted  by  the  CFD.  There  are  several  reasons  for  the  poor 
performaiKe  of  the  CFD.  and  which  continuing  research  is  attempting  to  resolve.  One  major  issue  is  the 
gnd  the  current  grid  is  loo  corrse  out  in  the  main  stream,  in  order  to  capture  wall  effects  (which  are  not 
being  studied  here).  This  has  been  identified  because  the  matched  pressure  case  was  fairly  well 
represented:  the  core  remained  next  to  the  wall  for  almost  the  entire  domain,  which  was  in  the  area  where 
the  gnd  density  was  highest.  The  overpressure  solution  core  moved  out  into  the  low  grid  density  region 

almost  immediately  The  grid  is  being  regenerated  to  have  better  main  flowlleld  resolution.  The  othc. 
adjustment  that  is  being  made  is  the  return  of  SCf  to  0.7,  where  binary  theory  indicates  it  should  be.  Also. 

grid  Mocking  issues  are  being  investigated. 

3.  CONCLUSIONS  AND  FUTURE  RESEARCH 

The  ability  of  CFD  solutions  to  accurately  model  complex  mixing  and  reacting  flowTiclds  is  the 
major  insure  that  has  been  addressed  here.  The  first  two  of  three  steps  have  been  undertaken.  The  first 
step  IS  to  determine  the  effectneness  of  GASP  in  modelling  reacting  flowficlds.  Here,  the  internal 
influences  of  modelling  were  examined.  The  effects  of  free  stream  turbulence  (and  also  different  models) 
was  very  dramatic  on  the  flame  speed  observed  in  a  shock-induced  combustion  problem.  The  effect  of 
turbulent  Schmidt  number,  which  is  the  parameter  which  directly  controls  the  amount  of  diflusion  in  the 
model,  was  examined  for  an  experimental  flowficld.  The  results,  in  terms  of  flame  location,  were 
improved  with  respect  to  the  experimental  results,  but  there  were  other  modelling  issues  that  need  to  be 
rcsoKed  for  that  problem  The  turbulent  Schmidt  number  does  have  a  large  impact  on  the  diffusion  in 
these  high  speed  reacting  cases,  so  that  the  models  can  be  tuned  for  specific  problems.  A  third 
comparison  was  done  between  different  chemistry  models  in  GASP,  for  a  shock-induced  combustion 
ptiblem  The  finite  rate  kinetics  turned  out  to  be  very  sensitive  to  the  model  chosen,  especially  for 
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problems  near  the  ignition  point.  For  a  general  combustion  problem.  GASP  modelled  the  expected 
physics  well. 

Mixing  in  supersonic  flowficlds  was  also  examined.  Two  1 5  degree  flush  wall  port  injectors  on  a 
flat  plate  were  modelled  and  the  mixing  results  were  compared  to  experimental  data  and  other 
computational  results.  The  ciurent  solution  over-diiluscd  the  core  and  overpr^icted  the  penetration  of 
the  core,  especially  for  the  overpressure  injection  case.  The  results  were  much  better  for  the  matched 
pressure  case.  This  effect  was  also  seen  in  other  solutions  for  this  problem  with  GASP.  The  SPARK  code 
generated  a  good  solution,  illustrating  that  the  problem  can  be  modelled  computationally.  The  major 
reason  for  the  discrepancy,  in  this  case,  is  postulated  to  be  the  grid.  Future  research  will  concentrate  on 
improving  the  mixing  solution. 

The  final  part  of  this  research  is  to  model  both  mixing  and  reacting  solutions  for  swept  ramp 
injectors  at  a  flight  Mach  number  of  13,  and  evaluate  the  solution  using  experimental  dau.  The  final  goal 
of  this  and  future  research,  is  to  demonstrate  GASP  as  a  tool  that  can  be  used  to  predict  the  performance 
of  sciamjet  combustors. 
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APPENDIX:  CHEMISTRY  MODEL  DETAILS 


1.  Dnimmond  Chemistry  Model:  7  species,  7  reactions. 
Species:  N2,  O2,  H2.  OH,  H2O.  O.  H 
Reactions: 

I.O2  +  H2  ==>OH-K)H 
3  H2  +  0H=>  H2O  +  H 
5.  OH  +  OH  ==>  H2O  +  O 
7.  H  +  H  +  N2  H2  +  N2 


2.  O2  +  H  =>  OH  +  O 
4.  H2  +  O  ==>  OH  +  H 
6.  OH  +  H  +  N2  =>  H2O  +  N2 
(N2  :  third  bod>'  for  ail  reactions) 


2.  Evans  and  Schexnayder  1  &1  Chemistry  Models: 

Model  1:  7  species,  8  reactions  (bold  face).  Model  2:  12  species.  25  reactions 
Species:  N2, 02,  H2,  OH,  H2O,  O,  H,  NO.  N02,  HO2,  HNO2,  N 
Reactions: 


1.  O2  +  N2 
3.  OH  +  N2 
5.  O2  +  H  • 
7.H2O  +  O 


2O  +  N2 
>  O  +  H  +  N2 
OH  +  O 
OH -I- OH 


2.  H2  +  N2  “■ 
4.  H2O  +  N2 
6.  H2  +  O  * 
8.  H2O  +  H 


->  2H  +  N2 
— >  OH  +  H  +  N2 
>  OH  +  H 
OH  +  H2 


9.  HNO2  +  N2  ==>  NO  +  OK  +  N2 
1 1.  HO2  +  N2  =>  H  +  O  +  N2 
I3.H2  +  O2  =*>H  +  H02 
I5.H2O  +  O  =>H  +  H02 
I7.H2O  +  O  =>0H  +  H02 
19.0  +  N  =>N  +  NO 
21.0  +  N0  ==>N  +  02 
23.NO  +  O2  =>0  +  N02 
25.  NO2  +  OH  NO  +  HO2 

3.  HASP  4  (Hydrogen/NOj  extension):  12  species,  24  reactions 

Species:  N2, 02,  H2,  NO,  OR  NO2,  HO2.  HNO,  H2O.  N.O,H 


10.  NO2  +  N2  =>N0  +  0  +  N2 
12  .  H2  +  O2  =>OH+OH 
14.  OH  +OH  ==>  H  +  HO2 
16.  OH  +  O  =>  O  +  HO2 
18.  H2O  +  OH  =>  H  +  HO2 
20.H  +  N  =>N  +  OH 
22.  NO  +  OH  ==>H  +  N02 
24.  NO2  +  H  =>  H  +  HNO2 


Reactions: 

I. OH  +  H2  =“>H  +  H20 
3.  O  +  H2  ==>  H  +  OH 
5.H+HO2  ==>OH  +  OH 
7.OH  +  HO2  ==>H20  +  02 
9.H  +  OH  +  N2  =>H20  +  N2 

II. H  +  O  +  N2  ==>oh  +  n 
13.  OH  +  OH  =>0  +  H20 
13.0  +  N2  =>N  +  N0 

I7.H  +  NO  +  N2  ==>HN0  +  N2 
19.  H  +  HNO  =>H2  +  N0 
21.  OH  +  HNO  ==>  H2O  +  NO 
23.  O  +  NO2  =>02  + NO 


2.H  +  O2  =>0  +  0H 
4.  H  +  HO2  =>H2  +  02 
6.  O  4  HO2  =>  OH  +  O2 
8.  H  +  O2  ==>  HO2 
10.  H  +  H  +  N2  ==>  H2  +  N2 
I2.O  +  O  +  N2  =>02  +  N2 
14.  O  +  NO  =>  N  +  O2 
16.  H  +  NO  =>  N  +  OH 
18.  O  +  NO  +  N2  =>  NO2  +  N2 
20. 0  +  HNO  =>  OH  +  NO 
22.  H  +  NO2  =>  OH  +  NO 
24.  HO2  +  NO  =>  OH  +  NO2 
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-  k-e  rtirbnlriirr  (Chirn']  MoJrl), 

Frfes»re«m  Turbuifnrr  intrtxiiy  •=  2.0*'. 


" -  k-£  Turbuifnce  (C'hifn'i  Model). 

Freestrcam  Turbulence  intensiiy  =  0.2*/. 


Figure  1 ;  Sketch  showing  the  ilame  fronts  computed  for  different  tuibulcnce  models 

and  model  geonMtry.  The  flame  front  is  taken  to  be  the  location  of  maximum 
water  mass  fraction  change. 


Figure  2;  Burrows  and  Kurkov  test  case,  computational  and  experimental  results. 

a),  geometry  sketch,  b).  inflow  conditions,  c).  current  results  for  water  mole 
fractions  compared  with  data  and  a  prexiously  published  CFD  solution.  Inset 
is  the  water  mole  fraction  Figure  published  by  Burrows  and  Kurkov.  showing 
that  the  original  analytical  solution  is  as  inaccurate  as  the  current  computation 
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Figure  5:  Quantitative  results  comparisons  between  the  chen)istr>  models. 


Figure  7;  Computational  results  for  hydrogen-air  reacting  flowfield.  Premixed  hydrogen- 
air  inflow  at  left,  over  a  backstop,  where  it  is  ignited  by  the  hydrogen-air  jet 
illustrated  in  Figure  6.  Upper  boundary  after  backstep  is  unconfined. 
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a).  Matched  pressure  injection 


Figure  10:  Matched  pressure,  X/D  =  80, 
helium  mass  ftaaion  contours. 


Figure  1 1 :  5X  overpressure,  X/D  =  80. 
helium  mass  fraction  contours. 


a)  SPARK  CFD  soluUon/experiment  b)  previously  published  GASP  solution 

comparison  comparison  with  experiment 


Figure  12:  5X  overpressure,  X/D  *  80.  helium  mass  fraction  contour  comparison 
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MULTIDIMENSIONAL  CONJUGATE  HEAT  TRANSFER  ANALYSIS 
FOR  THE  ARNOLD  ENGINEERING  DEVELOPMENT  CENTER 
HEAT-Hl  TEST  UNIT  NOZZLE 


Michael  A.  Weaver,  MSAE 
Graduate  Student 
School  of  Aerospace  Engineering 
Georgia  Institute  of  Technology 

Abatcact 

A  method  for  unsteady,  axisymmetric,  conjugate  heat  transfer  analysis 
was  developed.  The  conjugate  heat  transfer  domain  comprises  coflowing  high 
temperature  air  and  subcooled  water  coolant  on  opposite  sides  of  a  copper- 
zirconium,  converging  nozzle.  Heat  transfer  through  the  nozzle  wall  is 
characterized  by  solid  body  conduction  with  convection  boundary  conditions 
along  the  air  side  and  water  side  of  the  nozzle  wall.  The  air  side  heat 
transfer  is  characterized  by  forced  convection  with  a  turbulent  boundary 
layer.  The  water  side  heat  transfer  is  characterized  by  forced  convection, 
subcooled,  nucleate  boiling.  Convective  heat  transfer  coefficients  on  each 
side  of  the  nozzle  wall  are  functions  of  the  wall  temperature  and  the 
respective  flow  properties,  thus  coupling  the  three  regions  of  the  domain. 
The  solution  method  marches  in  time,  solving  at  each  time  step  for  the  nozzle 
wall  temperature  distribution,  the  flow  properties  on  each  side  of  the  nozzle 
wall,  and  for  the  convective  heat  transfer  coefficients.  The  algorithm 
terminates  when  either  the  steady  state  is  achieved  or  nozzle  wall  failure 
conditions  are  reached.  Preliminary  results  are  shown  for  run  conditions  at 
which  nozzle  wall  survival  has  been  experimentally  verified. 
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MULTIDIMENSIONAL  CONJUGATE  HEAT  TRANSFER  ANALYSIS 
FOR  THE  ARNOLD  ENGINEERING  DEVELOPMENT  CENTER 
HEAT-Hl  TEST  UNIT  NOZZLE 

Michael  A.  Weaver 

Introduction 

Within  the  context  of  this  study,  conjugate  heat  transfer  describes  the 
coupling  between  fluid  flow  over  a  structure  and  conductive  heat  transfer 
within  the  structure.  Such  coupling  requires  simultaneous  consideration  of 
both  the  convection  and  conduction  phenomena.  Conjugate  heat  transfer  occurs, 
for  example,  in  an  actively  cooled,  high  temperature  wind  tunnel  nozzle.  In 
this  case,  convective  heat  transfer  through  the  working  fluid  boundary  layer 
to  the  nozzle  wall  is  coupled  with  the  nozzle  wall  temperature.  Convective 
heat  transfer  from  the  nozzle  wall  to  the  coolant  is  also  coupled  with  the 
nozzle  wall  temperature.  In  turn,  the  conductive  heat  transfer  within  the 
solid  nozzle  material  is  coupled  through  convection  boundary  conditions  to 
flow  of  the  working  fluid  and  coolant  over  the  nozzle  wall. 

This  study  addresses  conjugate  heat  transfer  in  the  Arnold  Engineering 
Development  Center  (AEDC)  HEAT-Hl  Test  Unit  nozzle.  The  HEAT-Hl  Test  Unit 
(hereafter  referred  to  as  HEAT-Hl)  is  an  arc-heated,  free  jet  test  facility, 
providing  extremely  high  enthalpy  air  flow.  Flows  with  enthalpies  ranging 
from  2,000  to  8,500  Btu/lbm  and  pressures  ranging  from  20  to  115  atm  are 
routinely  produced  [1].  Mach  nuinbers  ranging  from  1.8  to  3.5  are  achieved 
with  interchangeable  nozzles.  Coflowing  water  with  sub-cooled,  nucleate 
boiling  provides  active,  backside  cooling  for  the  nozzle  wall. 

Proposed  future  applications  for  HEAT-Hl  require  pressures  up  to  200 
atm.  Unfortunately,  nozzle  wail  structural  failure  due  to  heat  load  occurs  at 
pressures  in  the  range  120  to  130  atm.  The  need  exists  to  predict  such  heat 
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This  and  the  previcua  studios  iacus  tho  HEAT-Hi  Mjch  1.8  r.oztie. 
This  axisintjnetric  nozzle  is  2.3  in,  long,  with  a  threat  aiaiTLeter  or  0.9  in. 
The  material  compositicn  is  copper-zirconium.  The  r.ezzie  ar.d  r:urtoundin.‘i 
cooling  jacket  are  approximated  in  Fig.  2.  Water  flows  Letwee.n  the  nozzle  and 
c  iing  jacket,  providing  backside  cooling  for  the  nozzle  wall.  2\way  from  the 
end  regions,  heat  transfer  througn  the  nozzle  wail  is  approximately  one- 
dimensional.  Near  the  flanged  ends  of  the  nozzle,  "he  solid  body  heat 
transfer  becomes  axisymraet ric  in  nature. 


Fig.  2.  Idealization  of  the  HEAT-Hl  nozzle  and  cooling  jacket. 

flow  over  the  nozzle  wall  is  assumed  steady,  compressible,  and 
turbulent,  with  a  high  temperature  boundary  layer.  Outside  the  boundary 
layer,  flow  is  assumed  to  be  steady,  compressible,  nonuniform,  dissociated, 
equilibrium  air,  with  decreasing  stagnation  temperature  approaching  the  wall. 
Transport  properties  for  the  high  temperature,  high  pressure  air  are  strong 
functions  of  temperature  and  weak  functions  of  pressure. 

Water  flow  between  the  coolrng  jacket  and  nozzle  wall  is  assumed  steady, 
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incorrpressibie,  viscous,  a:;a  turDuioi.c  .  Tbia  Aitot  -riCOio  cos  C:>:irs  f  r.cj  active 
cooling  region  in  the  subcocied  state.  Subcooled  rcuclearo  bcil:ng  ir  assuir.vd 
to  occur  at  the  nozzle  wall,  while  the  outer  cooling  jacket  wall  is  assurred 
adiabatic.  Transport  properties  for  the  subcoole-d  water  are  lunctions  of 

temperature  a.nd  pressure. 

Heat  load  structural  deformation  and  failure  in  the  HEAT-Hl  nozzles  has 
been  observed  in  the  region  between  the  nozzle  inlet  and  the  nozzle  throat. 
Nozzle  wall  failure  has  not  been  observed  downstream  cf  the  nozzle  throat. 
For  this  reason,  the  current  analysis  is  limited  to  the  HEAT-Hl  Mach  1.8 

nozzle  geom.etry  from  the  nozzle  inlet  to  the  nozzle  throat.  This 

simplification  leads  to  a  nonphysical  longitudinal  boundary  at  the  nozzle 
throat.  Due  to  the  approxim.ately  one-dimensional  (radial)  nature  of  heat 
transfer  at  the  throat,  this  longitudinal  throat  boundary  is  assum.ed 
adiabatic.  Tie  three  regions  of  the  conjugate  heat  transfer  domain  (nozzle, 
air,  and  water)  are  summarized  in  Fig.  3  with  the  adiabatic  boundary 

conditions  indicated. 


ADIABATIC  ////// 


Fig.  3.  Conjugate  heat  transfer  domain  for  the  HEAT-Hl  nozzle. 
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Matlaad-g^  Anaiyaia 


At  least  three  modes  of  heat  load  failure  exist  for  the  nozzle  wall. 
First,  the  temperature  in  the  nozzle  wall  may  reach  the  melting  point  of  the 
material,  causing  immediate  catastrophic  failure.  Second,  the  critical  heat 
flux  from  the  nozzle  wall  to  the  boiling  coolant  may  be  reached,  leading  to 
wall  "burnout".  Third,  the  time  span  and  magnitude  of  the  heat  load  may 
permit  plastic  deformation  great  enough  to  produce  structural  failure,  even 
though  wall  temperatures  remain  below  the  melting  point  of  the  material. 

The  failure  mode  being  examined  determines  the  appropriate  analysis 
method.  The  first  two  failure  modes  can  be  analyzed  with  either  a  steady  or 
an  unsteady  approach.  The  steady  conjugate  heat  transfer  problem  could  be 
solved,  and  then  the  solution  checked  for  nozzle  wall  temperatures  exceeding 
the  material  melting  point,  or  for  the  critical  heat  flux  being  reached. 
Alternatively,  the  unsteady  conjugate  heat  transfer  problem  could  be  marched 
in  time  until  either  the  nozzle  wall  temperatures  exceed  the  material  melting 
point,  the  critical  heat  flux  is  reached,  or  the  steady  state  is  obtained. 
The  third  failure  mode,  by  its  unsteady  nature,  must  be  analyzed  with  an 
unsteady  approach.  The  unsteady  method  already  described  could  be  used  for 
the  third  failure  mode,  with  the  inclusion  of  an  additional  step  to  determine 
structural  deformation,  and  a  check  for  structural  yielding. 

In  the  current  study,  only  the  first  two  failure  modes  are  considered, 
but  the  unsteady  approach  has  been  adopted  for  its  future  applicability  to 
plasticity  analysis  for  the  third  failure  mode.  The  algorithm  developed  is  as 
follows : 

a)  Assume  an  initially  constant  nozzle  wall  temperature  distribu¬ 
tion  equal  to  the  water  inlet  temperature. 
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b>  Calculate  the  initial  water  flow  properties  along  the  water 
side  and  the  initial  air  flow  properties  along  the  air  side  of 
the  nozzle. 

c)  Calculate  the  initial  heat  transfer  coefficients  for  the  water 
side  and  the  air  side  of  the  nozzle. 

d)  Use  the  heat  transfer  coefficients  and  effective  fluid  tempera¬ 
tures  for  the  air  side  and  water  side  of  the  nozzle  in  the  ini¬ 
tial  time  step  of  an  unsteady  heat  transfer  analysis  of  the 
copper-zirconium  nozzle  to  obtain  the  new  temperature  distribu¬ 
tion  . 

e)  Calculate  the  new  water  flow  properties  along  the  water  side 
and  the  new  air  flow  properties  along  the  air  side  of  the 
nozzle . 

f)  Calculate  the  new  heat  transfer  coefficients  for  the  water  side 
and  the  air  side  of  the  nozzle. 

g)  Use  the  heat  transfer  coefficients  and  effective  fluid  tempera¬ 
tures  for  the  air  side  and  water  side  of  the  nozzle  in  the  next 
time  step  of  an  unsteady  heat  transfer  analysis  of  the  copper- 
zirconium  nozzle  to  obtain  the  new  temperature  distribution. 

h)  Return  to  step  e) ,  until  either  the  steady  state  is  reached, 
the  critical  heat  flux  is  reached,  or  the  nozzle  wall  tempera¬ 
ture  exceeds  the  copper-zirconium  melting  point. 

At  each  time  step,  this  algorithm  calculates  the  steady  state  flow  and  heat 
transfer  properties  of  air  and  water,  assuming  the  wall  temperature 
distribution  is  in  thermal  equilibrium  after  each  time  step.  For  this 
assumption  to  remain  valid,  the  time  step  size  must  be  less  than  or  equal  to 
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the  characteristic  time  for  heat  transfer  in  the  nozzle  wail  material. 

This  approach  requires  numerical  models  (based  on  either  analytic 
methods  or  engineering  mathematical  correlations)  for  the  unsteady  solid  body 
heat  transfer  in  the  copper-zirconium  nozzle,  the  flow  properties  and  heat 
transfer  coefficients  on  the  air  side  of  the  nozzle,  and  the  flow  properties 
and  heat  transfer  coefficients  on  the  water  side  of  the  nozzle.  Each  of  these 
models  will  be  described  for  the  three  regions  of  the  conjugate  heat  transfer 
domain . 

(i)  Copper-zirconium  nozzle  wall; 

Unsteady,  axisymmetric,  solid  body  heat  conduction,  with  no  internal 
heat  generation,  and  with  temperature  dependent  thermal  conductivity  is 
governed  by 


aT  1  a  r.  aT^  a  f,  a? 


pCo  —  =  —  —  kr  —  +  —  kr  — — 
^  dt  r  ar  V  ar  J  ay  dy 


The  adiabatic  wall  boundary  condition  is  given  by 


The  convection  boundary  condition  is  given  by 


Kill 


Here,  =  isobaric  specific  heat, 


=  thermal  conductivity. 


=  heat  transfer  coefficient, 

=  boundary  normal  coordinate. 


=  radial  coordinate. 


=  t ime , 


solid  body  temperature. 
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eCt 


=  effective  fluid  temperature, 


=  wall  temperature, 
y  =  longitudinal  coordinate, 

p  =  mass  density. 

The  finite  element  program  TR.^X  [3]  is  used  to  solve  this  problem.  At 
each  time  step,  the  heat  transfer  coefficients  and  effective  fluid 
temperatures  are  specified.  The  effective  fluid  temperature  for  the  water 
side  is  the  local  static  temperature.  For  the  air  side,  the  effective  fluid 
temperature  is  the  local  adiabatic  wall  temperature.  Using  these  values, 
program  TRAX  reads  the  nodal  temperature  distribution  from  the  previous  time 
step  and  calculates  the  new  nodal  temperature  distribution.  The  finite 
element  mesh  in  Fig.  4  shows  the  geometry  for  the  HEAT-Hl  Mach  1.8  nozzle. 


Fig.  4.  Finite  element  model  of  the  HEAT-Hl  Mach  1.8  nozzle. 

Material  properties  for  elemental  copper  were  used  for  the  copper- 
zirconium  nozzle.  Values  for  Cp,  k,  and  p  were  linearly  interpolated  between 
values  specified  at  491.67°  r  and  2,500°  R.  The  melting  point  for  elemental 
copper  is  approximately  2,410°  R. 
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(11)  Air  side  of  the  nozzle  wall: 

Steady  state  heat  transfer  through  a  turbulent,  compressible  boundary 
layer  in  accelerating,  axisymmetric  duct  flow  is  given  by  the  method  of  Ambrok 
[4], 

h 

GCp  ■ 

Here,  c^  =  isobaric  specific  heat, 

G  =  mass  flux, 

h  =  heat  transfer  coefficient, 

Pr  =  Prandtl  number, 

R  =  wall  x.adius, 

St  =  Stanton  number, 

=  adiabatic  wall  temperature  =  +  Pr'''’(T^  -  T^) , 

=  stagnation  temperature, 

T,  =  static  temperature, 

=  wall  temperature, 

X  =  wall  arc  length  at  position  of  interest, 

=  absolute  viscosity, 

4  =  variable  of  integration  (wall  arc  length)  . 

This  equation  is  based  on  solving  the  energy  integral  equation  for  the 
boundary  layer  [5]  .  It  is  valid  for  smoothly  varying  wall  temperature  and 
assumes  a  cooled  wall.  The  boundary  layer  is  assumed  to  originate  at  the 
lower  limit  of  the  integral.  All  properties  are  evaluated  at  the  local  static 
temperature . 

The  flow  properties  for  dissociated,  equilibrium  air,  as  required  for 
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the  HEAT-Hl  application,  are  obtained  by  mathematical  correlation.  An 
effective  gas  constant  and  specific  heat  ratio  are  used  in  the  quasi-one- 
dimensional,  isentropic  flow  equations  to  obtain  mass  flux,  static 
temperature,  and  isobaric  specific  heat.  The  effective  constants  are  chosen 
such  that  the  correlated  flow  variables  compare  well  with  flow  variables  from 
a  predetermined  equilibrium  gas  solution.  The  boundary  layer  edge  stagnation 
temperature  from  the  real  gas  solution  is  used  as  the  effective  stagnation 
temperature  in  the  isentropic  flow  equations.  This  value  differs  from  the 
nominal  HEAT-Hl  run  condition  stagnation  temperature  due  to  nonuniformity  of 
the  flow.  The  nominal  HEAT-Hl  run  condition  stagnation  pressure  is  used  as 
the  effective  stagnation  pressure  in  the  isentropic  flow  equations. 

The  air  transport  properties  are  obtained  from  equilibrium  gas  tables 
[6]  .  Absolute  viscosity  and  Prandtl  number  are  bilinearly  interpolated  from 
values  tabulated  as  functions  of  pressure  and  temperature.  The  local  pressure 
and  temperature  used  for  interpolation  are  obtained  from  the  flow  properties 
correlation. 

The  heat  transfer  coefficients,  h,  and  adiabatic  wall  temperatures, 
are  determined  along  the  air  side  of  vhe  nozzle  wall  using  these  models. 
These  values  are  then  used  to  specify  the  air  side  boundary  conditions  for  one 
time  step  of  the  finite  element  analysis. 

(ill)  Water  side  of  the  nozzle  wall: 

Steady  state  heat  transfer  during  forced  convection,  subcooled  nucleate 
boiling  in  an  annular  passage  is  predicted  by  the  correlation  of  Shah  [7,8,9]. 
This  correlation  assumes  the  total  convective  heat  flux  is  the  sum  of  the 
single-phase  convection  heat  flux  and  the  nucleate  boiling  heat  flux, 

q  =  q  +  q  , 

*  “spc 
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where,  q 


=  two-phase  convection  heat  flux 


=  single-phase  convection  heat  flux, 
q^j,  =  nucleate  boiling  convection  heat  flux. 

The  single-phase  heat  flux,  for  turbulent  flow,  is  determined  from  the 


where,  =  equivalent  annulus  diameter, 

G  =  mass  flux. 


h^p  =  single-phase  convection  heat  transfer  coefficient 

k  =  thermal  conductivity, 

Pr  =  Prandtl  number, 

T,  =  static  temperature, 

=  wall  temperature, 

H.  =  absolute  viscosity. 

All  quantities  are  evaluated  at  the  local  static  temperature. 

The  nucleate  boiling  convection  heat  flux  is  determined  from  the 


correlation  of  experimental  data  for  fully  developed  flow  boiling.  This 
follows  the  recommendation  of  Bergles  and  Rohsenow  [10]  when  they  demonstrated 
the  nucleate  boiling  convection  heat  flux  should  not  be  determined  with  a  pool 
boiling  correlation. 

The  Shah  correlation  for  forced  convection,  subcooled  nucleate  boiling 
in  an  annular  passage  is  expressed  as, 

q  '  -  T,,JMax(230Bo^^^  1)  +  a(T„^  -  T,)  ]  =  h,p(T„^^^  -  T,)  , 

where.  Bo  =  boiling  number  =  q/(i,^G), 
ij^  =  heat  of  vaporization. 
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h.p  =  two-phase  convection  heat  transfer  coefficient, 

=  saturation  temperature  at  local  static  pressure. 

Also, 

fO,for  fully  developed  boiling 
=  1 

[l,for  local  or  partial  boiling 
Here,  local  or  partial  boiling  occurs  when, 

>2  or  >  6.  3xl0‘  Bo'^' , 


(t,,,  -  tJ 


Otherwise,  fully  developed  boiling  is  assumed. 

If  no  boiling  is  present  {i.e.  ,  then  the  heat  transfer 

coefficient  is  chat  for  single-phase  convection,  h  =  h^p.  Otherwise,  the  heat 
transfer  coefficient  is  that  for  two-phase  convection  from  the  Shah 
correlation,  h  =  h^^,  which  already  includes  a  contribution  from  single-phase 
convection . 

The  nucleate  boiling  convection  heat  flux  predicted  by  the  Shah 
correlation  {and  subsequently  used  to  obtain  the  two-phase  convection  heat 
transfer  coefficient)  must  be  compared  to  a  predicted  critical  heat  flux.  The 
critical  heat  flux  model  used  is  derived  from  the  Rousar-Chen  model  [11]  for 
flat  plates.  Fred  Shope  at  Calspan  modified  the  Rousar-Chen  critical  heat 
flux  by  correlating  data  from  curved  plate  flows  to  account  for  transverse 
flow  acceleration  effects  due  to  the  curvature.  The  Shope-Rousar-Chen 
correlation  for  critical  heat  flux  in  curved  plate  flows  is  given  by 

'?chi  =  [^30  +  131.  778(l  -  e'  -  T,  )]j^l  +  0.  2(1  -  e"  )  , 

where,  g  =  gravitational  acceleration  (ft  s‘^), 

=  critical  heat  flux  (Btu  ff^  s'=), 

R  =  radius  of  curvature  (ff'). 
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=  Static  temperature  (“R)  , 

=  saturation  temperature  at  local  static  pressure  (°R)  , 

V  =  flow  velocity  (ft  s"')  . 

The  flow  properties  of  water  through  the  annular  passage  are  determined 
by  application  of  the  steady  state  conservation  laws  to  a  one-dimensional 
control  volume  with  area  change.  A  control  volume  is  assumed  with  the  known 
quantities, 

=  inlet  area, 

A^^^  =  exit  area, 

A_^^  =  wetted  surface  area, 

=  equivalent  annulus  diameter, 
fa  =  mass  flow  rate, 

pj„  =  inlet  static  pressure, 

Q  =  heat  transfer  rate  into  control  volume, 

=>  inlet  static  temperature, 

^  =  absolute  viscosity, 

p  =  constant  mass  density. 

Now,  from  conservation  of  mass  in  the  control  volume,  the  inlet  and  exit 
velocities  must  be 


V,„  = 


m 


,  and  = 


m 


A,,,p 

The  expression  for  conservation  of  momentum  in  the  control  volume  is 
given  by 


-  PoutAcut  +  (—  ](Aou.  -  A.n)  -  -  V,  J . 

Here  the  turbulent  wall  shear  stress  t  can  be  determined  with  Prandtl's 
universal  law  of  friction  for  smooth  pipes  [12],  and  the  known  values 
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4,  and  p.  This  leaves  exit  pressure  as  the  only  unknown  in  the  momentum 
equation . 

Conservation  of  energy  for  the  control  volume,  assuming  no  work  being 
done,  is  given  by 


where,  h,^  =  inlet  specific  enthalpy, 

h^^j^  =  exit  specific  enthalpy. 

This  equation  is  used  in  combination  with  the  known  values  p.„,  p,^^,  T.^,  and 
the  thermodynamic  relations  for  water, 

to  solve  for  the  exit  static  temperature  of  the  control  volume. 

Note  that  the  Shah  correlation  and  the  method  for  determining  local  flow 
properties  are  coupled,  due  to  the  equivalence  of  the  two-phase  convection 
heat  flux  and  the  control  volume  heat  transfer  rate  per  unit  heated  area.  For 
this  reason,  iteration  between  determining  the  local  flow  properties  and 
determining  the  two-phase  convection  heat  flux  is  required. 

The  water  transport  properties  are  obtained  from  the  standard 
thermodynamic  relations  for  subcooled  water  [13,14]. 

The  heat  transfer  coefficients,  h,  and  local  static  temperatures,  T^, 
are  determined  along  the  water  side  of  the  nozzle  wall  using  these  models. 
These  values  are  then  used  to  specify  the  water  side  boundary  conditions  for 
one  time  step  of  the  finite  element  analysis. 

Beauita 

Using  this  method  of  analysis,  an  axisymmetric  conjugate  heat  transfer 
solution  for  the  HEAT-Hl  Mach  1.8  nozzle  was  obtained.  The  Mach  1.8  nozzle 
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has  been  shown  to  survive  at  the  run  conditions  used  in  this  analysis.  These 
conditions  are. 

Air  total  pressure;  126.5  atm 

Air  total  temperature:  9,000°  R 

Air  total  enthalpy:  3,480  Btu/lb, 

Water  mass  flow  rate;  11.54  Ib./s 

Water  inlet  temperature;  557°  R 

Water  inlet  pressure:  1,000  psia 

The  effective  properties  for  air  were  determined  to  be 
Gas  constant;  0.07400  Btu  Ib^'-  °R'* 

Specific  heat  ratio;  1.24 

Total  temperature:  7,513°  R 

Also,  a  constant  Prandtl  number  of  0.72  was  used  for  the  air  side.  The 
effective  total  temperature  corresponds  to  the  boundary  layer  edge  total 
temperature  for  the  given  run  conditions  in  the  Mach  1.8  nozzle.  This 
effective  total  temperature  (not  the  nominal  total  temperature  of  9,000°  R) , 
was  used  for  the  determination  of  flow  properties  on  the  air  side  of  the 
nozzle . 

The  algorithm  was  marched  in  time  with  a  step  .size  of  0,00001  seconds. 
After  a  solution  time  of  0.065  seconds  (6,500  time  steps),  the  1-norm  of  the 
nodal  temperature  change  dropped  from  4.2%  to  0.18%.  The  maximum  nodal 
temperature  change  dropped  from  0.54%  to  0.0059%.  At  this  point,  neither  the 
nozzle  wall  melting  temperature,  nor  the  criuical  heat  flux  along  the  water 
side  of  the  nozzle  had  been  reached. 

Results  from  this  solution  are  compared  to  results  from  the  Calspan  one¬ 
dimensional  analysis.  Nozzle  wall  temperature  distributions  for  the  air  side 
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included  in  the  one-dimensional  analysis.  The  .t.dsc  striking  differences 
between  the  one-dimensional  analysis  and  r.ho  axisym.T.et  ric  analysis  occur  near 
this  region  of  the  nozzle.  The  abrupt  temperature  change  of  the  a-Kisymmetric 
result  in  Fig.  6  demonstrates  t'le  influence  of  the  multidimensional  thermal 
effects  which  cannot  be  addressed  by  a  one-dimensional  approach.  This 
rapidity  of  temperature  change  also  indicates  that  the  admittedly  coarse 
finite  element  mesh  should  be  refined  in  this  region. 

Conclusion 

The  development  of  a  method  for  unsteady,  axisyirmet ric,  conjugate  heat 
transfer  analysis  has  been  accomplished,  and  preliminary  results  indicate  its 
feasibility.  However,  useful  application  of  this  method  requires  refinement 
of  its  individual  components,  and  further  investigation  of  experimentally 
verifiable  test  cases.  The  space  restriction  of  this  report  prevents 
mentioning  more  than  a  few  potential  improvements. 

The  current  flow  model  for  the  water  side  of  the  nozzle  should  be 
replaced  with  an  axisymmetric,  incompressible,  viscous  flow  solver.  This  is 
required  to  ascertain  the  degree  of  recirculation  as  the  flow  winds  between 
the  nozzle  flange  and  cooling  jacket  at  the  water  inlet.  Refer  back  to  Fig.  3 
for  the  flow  path  of  water  in  the  HEAT-Hl  Mach  1.8  nozzle. 

The  effects  of  finite  element  mesh  refinement  and  solution  time  step 
size  must  be  investigated.  At  present,  only  one  finite  element  mesh  (see  Fig. 
4)  has  been  used.  This  mesh  probably  represents  an  upper  limit  on  coarseness. 
Conjugate  heat  transfer  solutions  should  be  obtained  with  progressively  finer 
meshes  to  understand  the  convergence  behavior  of  the  solution  method.  Also,  a 
characteristic  time  for  heat  transfer  in  the  nozzle  geometry  has  not  been 
adequately  investigated. 
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SETTING  CRITERIA  FOR  HALON  REPLACEMENT  AGENTS 


Timothy  Keen 
Graduate  Student 
Fire  Research  and  Testing  Center 
University  of  Florida 

Abstract 

The  program  to  replace  CFC^s  in  the  U.S.  Air  Force  inventory 
entails  the  examination  of  alternatives  for  the  various  major  uses: 
refrigerants,  solvents,  and  fire  suppressants.  Classic  decision 
matrices  allow  a  single  evaluator  to  both  rate  and  provide  weights 
for  each  criterion  against  the  various  alternatives.  The 

methodology  described  in  this  paper  allows  multiple  evaluators  to 
rank  the  criteria  in  order  to  generate  criteria  weights.  A 
symbolic  scheme  to  state  the  relative  importance  of  the  criteria 
and  a  system  for  "collapsing”  the  rankings  are  described. 
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SETTING  CRITERIA  FOR  HALON  REPLACEMENT  AGENTS 


Timothy  Keen 


INTRODUCTION 

The  Montreal  Protocol  of  1987  and  the  Clean  Air  Act  Amendments 
of  1990  have  dictated  that  the  class  of  chemicals  known  as 
chlorofluorocarbons  (CFC's)  be  banned  from  production  by  January  1, 
1995  and  from  use  by  January  1,  2000.  Largely  used  as 
refrigerants,  CFC's  also  play  and  important  role  as  firefighting 
agents  where  clean  fire  suppression  is  important.  In  computer  and 
communications  facilities  it  is  desirable  to  extinguish  the  various 
classes  of  fires  that  may  occur  with  minimal  disruption  to  ongoing 
operations  and  without  adding  to  the  damage  caused  by  the  fire 
suppression  method.  Use  of  water  sprinklers,  foams,  or  dry 
chemicals  necessitate  extensive  and  expensive  cleanup  operations 
after  a  fire.  A  subset  of  CFC's  known  as  Halons  has  provided  clean 
fire  suppression  capability  for  over  two  decades.  Halon  1301  has 
been  the  clean  agent  most  frequently  employed  to  protect 
computer/communications  facilities  from  fire  damage  and  collateral 
damage  that  would  be  a  function  of  the  fire  suppression  method.  A 
replacement  agent  for  Halon  1301  is  being  sought  to  met  the 
requirements  of  the  Clean  Air  Act  Amendments  of  1990  and  a  set  of 
criteria  and  a  decision  analysis  matrix  will  be  necessary  to  select 
the  follow-on  agent. 
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DISCUSSION 


In  another  U.S.  Air  Force  program  to  select  a  replacement 
agent  for  Halon  1211,  used  in  flightline  firefighting  and  in 
facility  and  aircraft  portable  firefighting  units,  a  decision 
analysis  matrix  was  utilized  rather  late  in  the  program  to  provide 
a  framework  to  justify  the  selection  of  the  replacement  agent. 
Halon  1211  was  compared  to  PFC-614  and  HCFC-123  to  determine  which 
of  the  replacement  agents  should  be  selected  for  further  testing. 
Criteria  such  as  agent  effectiveness,  acute  toxicity,  system 
conversion  costs,  purchase  cost,  ozone  depletion  potential  (ODP) , 
and  greenhouse  warming  potential  (GWP)  were  used  as  the  basis  for 
agent  evaluation.  The  decision  analysis  method  used  in  this 
application  had  several  undesirable  features  that  need  to  be 
corrected  for  future  programs  such  as  the  one  that  will  replace 
Halon  1301  as  the  total  flood  agent  for  occupied  facilities. 

First,  each  criterion  in  the  decision  matrix  was  treated 
identically  in  terms  of  weight.  Fire  suppression  efficiency 
received  the  same  weight  or  emphasis  as  toxicity  and  cost.  For 
all  practical  purposes  it  could  be  said  that  no  weighting  scheme 
was  utilized. 

Second,  the  scoring  of  each  alternative  was  accomplished  by 
awarding  three  points  for  the  alternatives  with  the  best 
performance  against  each  criterion  and  one  point  against  each 
criterion  showing  the  worst  performance.  This  spread  of  points  is 
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not  enough  to  significantly  differentiate  good  performance  from 
poor  perf oirmance . 

Third,  several  criteria  are  actually  counted  a  number  of 
times.  The  category  entitled  future  regulatory  phase-out  is 
accounted  for  several  times  since  OOP,  GWP,  HCFC's,  EPA  SNAP 
approval  are  all  connected  with  this  topic.  Consequently  out  of  16 
criteria,  four  are  counted  against  possible  phase-out. 

This  example  illustrates  several  of  the  pitfalls  associated 
with  decision  matrix  methods.  A  set  of  clearly  defined  criteria 
are  necessary  in  order  to  set  up  the  decision  matrix.  Each 
criterion  should  appear  once  in  the  matrix  and  there  should  be  no 
interdependence  of  criteria.  This  latter  recommendation  may  not 
always  be  possible  to  implement.  Nonetheless  it  should  be  utilized 
to  the  maximum  extent  possible.  A  weighting  scheme  should  also  be 
employed  in  order  to  indicate  the  relative  priority  of  the 
criteria.  The  weights  should  be  generated  using  input  from 
several  experts  in  the  field  in  order  to  minimize  the  influence  of 
any  single  rater  on  the  final  weights  generated. 

It  is  this  latter  point  that  is  perhaps  the  most  difficult  to 
achieve.  There  is  not  at  the  present  time  any  method  that  allows 
the  opinions  of  several  experts  to  be  synthesized  into  a  single 
outcome.  In  order  to  accomplish  this  synthesis,  a  method  for 
aggregating  the  inputs  of  multiple  experts  is  required. 
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METHODOLOGY 


The  approach  to  solving  the  problem  of  allowing  multiple 
experts  to  generate  criteria  weights  was  to  create  a  system  in 
which  experts  could  provide  their  inputs  in  a  symbolic  fashion. 
The  symbolic  system  consisted  of  listing  the  criteria  from  left  to 
right  with  the  most  important  criteria  being  in  the  leftmost 
position  and  the  least  important  in  the  far  right  position.  The 
following  symbols  are  used  to  define  the  relationship  of  the  left 
criterion  to  its  right  neighbor: 

=  the  criteria  are  of  about  equal  importance 
>  the  left  criterion  is  slightly  more  important 
>>  the  left  criterion  is  more  important 
»>  the  left  criterion  is  far  more  important 

For  a  scheme  with  12  criteria,  the  following  is  an  example  of 
how  a  single  evaluator  might  rank  the  criteria: 

5>6=4»7»>1=2>3»8»>9=10>11»13=12 

Note  that  each  criterion  is  used  only  once  and  that  all 
criteria  are  ranked. 


% 
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The  criteria  that  will  be  used  in  the  Halon  1301  program  are 
as  follows: 

*  Fire  Suppression  Efficiency 

*  Low  Residue  Level 

*  Low  Electrical  Conductivity 

*  Low  Metals  Corrosion 

*  High  Metals  Compatibility 

*  Stability  under  Long  Term  Storage 

*  Low  Toxicity 

*  OOP 

*  GWP 

*  Cost 

*  Production  Availability 

*  Extinguishment  Concentration 

*  Conversion  Cost  of  Facility 

In  order  to  determine  or  create  a  suitable  decision  analysis 
system,  these  criteria  must  be  analyzed  and  ranked  by  a  number  of 
experts.  The  rankings  created  by  the  "experts'*  will  then  be 
aggregated  into  a  single  expression  of  relative  importance  of  the 
criteria,  A  weight  will  be  generated  for  each  criterion  by  virtue 
of  the  relative  importance  of  the  criterion.  Each  alternative  will 
be  scored  on  a  relative  basis  against  each  criterion,  the  weights 
will  be  applied,  and  a  total  score  will  be  calculated.  Sensitivity 
analysis  will  be  utilized  to  determine  the  affects  of  the  weighting 
scheme  on  the  outcomes.  Finally  the  agents  with  the  greatest 
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number  of  points  will  be  selected  for  advanced  testing. 


CONCLUSIONS 

The  decision  analysis  matrix  system  described  in  this  report  is  a 
new  approach  to  generating  weights  for  a  wide  variety  of 
applications.  The  application  used  to  test  this  method  is  the 
selection  of  a  Halon  1301  replacement  agent.  However  the  main 
outcome  is  that  the  general  principles  described  herein  can  be  used 
for  any  case  that  would  benefit  from  the  inputs  of  multiple  experts 
in  the  generation  of  criteria  weights. 
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ABSTRACT 


The  effects  of  high  temperature  ‘- =  1000  jet  blast  on  runway 

surfaces  has  become  a  significant  concern,  particularly  with  the 
increasing  prominence  of  V/STOVL  ( vert ica I /shor t  take-off  and  landing) 
aircraft.  In  an  attempt  to  model  vertical  jet  impingement  heat 
transfer,  a  computer  code  has  been  developed  based  on  the  Hiemenz 
solution  of  the  Navier-Stokes  equations.  The  primary  use  of  the  code 
would  be  to  provide  values  of  heat  flux  and  wall  temperature  for  use 
as  input  to  finite  element  solid  mechanics  modeling  codes,  which  are 
currently  being  used  by  the  Jet  Blast  Research  Group  at  Tyndall  AFB 
Florida,  to  predict  stresses  in  pavement  materials  as  a  result  of 
high  heat  flux.  Preliminary  results  show  the  code  to  be  in  good 
agreement  with  experimental  data  and  analytical  calculations.  A 
research  project  is  being  proposed  which  would  continue  development 
of  the  computer  program. 
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3ACKGFQUND 

In  -"ecent  years,  significant  effort  has  oeen  airectea  tcwaras  fit-' 
aroblem  of  concrete  runway  deqraaaticn  as  a  result  o»  ceriodtc, 
intense,  heat  output  from  military  jet  engines.  wif^  the  inc'easinq 
prominence  of  V/STOVL  technology,  the  problems  associated  with  this 
phenomenon  nave  become  more  acute. 

As  with  most  research-or  lented  egineennq  endeavors,  the  solution  to 
this  problem  has  been  sought  using  both  eMperimental  and  analytical 
approaches.  Experiments  have  focused  on  testing  concrete  pavement 
?both  on-site  and  in  the  lab)  in  order  to  determine  the  mecharusms  of 
material  breakdown  as  a  result  of  thermal  input.  Full-scale  high  heat 
( appox  imatel  y  1000  "^F)  conditions  have  been  produced  using  actual 
aircraft  for  experimental  purposes. 

Analytical ly,  the  focus  has  been  on  modeling  the  mechanisms  which 
lead  to  concrete  failure  due  to  thermal  stresses,  with  the  more 
complex  models  considering  the  inhomogenous  nature  of  concrete  and  the 
cosequential  differential  thermal  stresses  created  by  heating.  Finte 
element  codes  are  used  mostly  in  these  efforts.  In  terms  of  mode' ing 
the  actual  thermal  (convective)  output  of  jet  engines,  and  the 
resultant  temperature  and  heat  flux  generated  when  directed  towards  a 
surface,  many  of  the  codes  of  this  type  (Bose  CM,  and  Aoeihoff  e£  at 
121,  for  example)  are  Nav ler-Stokes  solvers.  Nav ler-Stokes  codes, 
however,  are  inherently  complex  and  require  expensive  computer 
hardware  (supercomputers  mainly)  in  order  to  obtain  accurate  solutions 
in  a  short  time  period. 
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'HE  CODE 


Due  in  part  to  toe  cans iderat lonp  expressed  above,  a  arort.  easy 
to  use  jet  impingement  heat  transfer  code  has  bee  written  to 
predict  temperature  and  heat  flux  when  modeling  normal  jet  impingement 
on  a  solid  surface.  The  code  is  based  on  the  Hiemenz  stagnation 
region  solution  of  the  Navier-Stokes  equations  as  presented  by  White 
C33.  The  governing  differential  equation  is  solved  in  an  iterative 
manner  using  a  fourth-order  Runqe-Kutta  numerical  integration 
algorithm  coupled  with  an  interpolation  scheme  based  on  the 
half-interval  method.  Surface  temperatures  are  calculated  using  an 
equation  derived  by  Abelhoff  et  a.L  based  on  steady-state  conditions. 

The  primary  function  of  the  code  will  be  to  provide  heat  flux  and 
temperature  data  which  could  be  used  as  input  for  codes  which  model 
thermal-induced  stresses  in  concrete.  Another  possible  use  would  be 
as  a  means  of  calculating  the  minimum  thickness  of  runway  protective 
coatings.  At  this  time,  materials  are  being  tested  and  evaluated  for 
use  as  coatings  for  cocrete  pavement,  and  a  method  for  determining  the 
minimum,  or  critical  thickness  of  the  materialls)  would  be  beneficial 
from  an  economic  point  of  view. 
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Preliminary  results  oncairtea  frcm  the  coae  are  encGuraqino, 
that  agreement  has  been  achieved  with  both  experimental  '"^abie  1  aiia 
analytical  results  C4].  Particularly  significant  is  the  'act 
that  the  free  jet  itself  is  not  modeled;  Known  values  of  the  jet  e^it 
temperature,  exit  velocity,  along  with  the  nozzle  diameter  and  height 
are  all  that  are  required  as  input.  With  these  values,  the  jet 
velocity  and  temperature  near  a  surface  some  distance  from  the  nozzle 
can  be  calculated  using  the  general  equation  [53 


where  P  represents  the  desired  parameter  to  be  calculated,  C  is  c 
constant  being  equal  to  0.9  when  calculating  velocity,  and  0.B5  for 
calculating  temperature.  The  values  of  the  constant  were  derived 
using  results  from  calculations  based  on  the  method  of  [43  as  a  guide. 
The  variables  D  and  2  are  the  nozzle  velocity  and  height,  respectively. 

At  this  time,  the  code  predicts  steady-state  values  of  heat  flux 
and  surface  temperature.  Future  efforts  aimed  at  refining  the  code  will 
include  modifications  to  allow  calculation  of  transient  values  of  the 
parameters  mentioned  above.  Also,  an  attempt  will  be  made  to  model 
oblique-jet  heat  transfer  as  well.  Finally,  the  model  will  be 
extended  to  the  wall  jet  region  of  an  impinging  flow  in  order  to 
predict  heat  transfer  away  from  the  stagnation  point.  A  Research 
Initiation  Proposal  is  being  prepared  at  this  time,  and  subject  to 
approval  by  AFOSR,  these  objectives  will  he  the  focus  of  the  ''esearch 
effort . 
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SUMMARY 


A  computer  code  nas  been  developed  to  model  jet  impingement  heat 
transfer,  in  order  to  provide  thermal  input  conditions  for  finite 
element  codes  used  to  predict  the  effects  of  thermal  stresses  on 
runway  surfaces.  Initial  runs  have  resulted  in  data  which  is  in  good 
agreement  with  experimental  data  and  calculations  using  analytical 
methods.  Research  is  proposed  which  would  continue  development  and 
refinement  of  the  code. 
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TABLE  I 


Measured  versus  predicted  temperatures 
from  F-18  Auxilary  Power  Unit  impingement  tests  -Ref.  o) 


Mode 

Measured  Surface  Temp. 

Predicted 

T emp . 

Error 

MES  (Sch. 

I ) 

32E  ""f 

323.2 

"F 

0.4’/. 

ECS  (Sch. 

I ) 

328  "’f 

320.6 

r. 

F 

2 . 3'/. 

ECS  (Sch. 

II ) 

345  '^F 

320.6 

r 

7.17. 

* 

Note:  Data  from  production  nozzle  tests. 
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VORTICITY  AND  VELOCITY  MEASUREMENTS 

IN 

TRANSIENT  OSCILLATORY  SEPARATING  BOUNDARY  LAYER  FLOWS 


B.  Terry  Beck 
Associate  Professor 

PauiK..  Berg 
Graduate  Student 


Department  of  Mechanical  Engineering 
Kansas  State  University 


Abstract 

The  velocity  and  vorticity  distribution  within  a  transient  oscillatory  separating  boundary 
layer  was  investigated  using  a  single-component  Laser  Doppler  Velocimeter  System.  The  flow 
was  initiated  above  a  flat  plate  test  model  by  means  of  a  computer-controlled  rotating  spoiler 
(flap),  mooated  above  the  model  surface.  The  tests  were  conducted  in  a  water  tunnel  test 
facility,  and  dye  injection  was  also  utilized  for  visualization  of  the  flow  separation  phenomena. 
The  rotating  spoiler  subjected  the  plate  below  to  a  time-dependent  spatial  pressure  gradient, 
inducing  periodic  flow  separation  and  vortex  shedding  from  the  region  near  the  plate  and 
downstream  of  the  spoiler.  Measurements  of  both  horizontal  and  vertical  velocity  components 
were  made  by  rotating  the  optics  of  the  LDV  system.  These  profile  measurements  were  obtained 
for  discrete  angular  flap  positions,  thus  mapping  out  the  spatial  and  time-dependent  flow  field 
downstream  of  the  flap.  From  the  separate  velocity  component  profiles,  a  computerized  scanning 
algorithm  was  implemented  to  obtain  both  scan-averaged  velocity  and  velocity  gradient  fields. 
Using  this  technique  resulted  in  remarkably  smooth  results,  in  spite  of  the  limited  spatial 
resolution  of  the  transient  measurements.  Clear  evidence  of  reverse  flows  and  flow  bifurcation  is 
indicated  from  the  measurements  near  the  region  of  boundary  layer  separation.  The  effect  of  flap 
firequency  on  the  separation  phenomena  was  also  investigated. 
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MULTIVARIABLE  TRANSFER  FUNCTIONS  AND  OPTLMAL  PASSIVE 
DAMPING  FOR  A  SLEWING  PIEZOELECTRIC  LAMINATE  BEAM  ' 

Thomas  E.  Alberts  k  Travis  V.  DuBois 
Department  of  Mechanical  Engineering  and  Mechanics 
Old  Dominion  University 
Norfolk,  Virginia  23529-0247,  USA 
Phone;  (804)  683-3736 
Email:  taiberts(§lmem. odu.edu 

Abstract 

This  report  presents  the  development  and  experimental  verification  of  a  distributed  parameter 
model  fcr  a  slewing  beam  system  with  piezoelectric  actuators  and  sensors.  The  beam  is  pinned 
at  the  proximal  end,  an  endpoint  motion  sensor  is  attached  at  the  distal  end.  and  patches  of  thin 
piezoelectric  laminates  attached  to  its  surface.  The  differentia!  equation  of  motion  for  this  system 
is  transformed  to  Laplace  domain  treinsfer  functions  after  application  of  the  appropriate  boundary 
conditions.  Transfer  functions  relating  the  various  actuator /sensor  pairs  are  developed.  The 
transfer  functions  are  rationalized  using  a  Maclaurin  series  expansion  so  that  there  is  no  need  to 
assume  mode  shapes.  Experimental  results,  which  verify  the  model,  arc  presented  using  a  beam 
experiment  at  the  US  Air  Force  Academy,  Frank  J.  Seiler  Research  Laboratory.  The  existing 
clamped  beam  experiment  was  modified  through  the  addition  of  a  hinged  joint  and  appropriate 
instrumentation  to  carry  out  this  work. 

The  transfer  functions  are  eventually  to  be  used  to  develop  and  experimentally  validate  a 
simultaneously  optimal  active  and  passive  damping  design  for  the  experimental  system.  A  pre¬ 
liminary  damping  design  is  discussed  and  initial  experimental  results  presented. 


'  This  work  performed  in  coUaboraiion  with  Dr.  H.R.  Pota  of  the  Australian  Defence  Force  Academy. 
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A  NEURAL  NETWORK  MODEL 

OF  THE  UNSTEADY  AERODYNAMICS  ON  A  PITCHING  WING 

Willuun  E.  Fuller 
Research  Associate 
BioSer%'e  Space  Technologies 
University  of  Colorado,  Boulder 

ABSTRACT 

A  straight  wing  having  a  NACA  0015  cross-section  and  rectangular  phinform  was 
attached  to  a  circular  splitter  plate.  Starting  at  0  degrees  this  configuration  was  pitched  to  an 
angle  of  60  degrees  which  exceeded  the  static  stall  angle.  During  the  pitch  history  surface 
pressure  readings  were  obtained  from  15  pressure  transducers  spaced  between  0  and  90% 
chord.  A  total  of  54  data  records  were  obtained  which  covered  6  non-dimensional  pitch 
rates  (a"*")  ranging  between  0.0001  and  0.2  and  9  span  locations  ranging  between  0%  and 
80%  span.  These  unsteady,  vortex  dominated  flows  were  used  to  develop  an  artificial 
neural  network  (ANN)  model  of  the  unsteady  flow  field.  The  ANN  model  was  then  used  to 
mathematically  quantify  the  three-dimensional,  vortex  dominated,  unsteady  aerodynamics 
of  the  phenomenon.  A  linear  equation  system  (LES)  was  derived  from  the  weight  matrices 
of  the  ANN.  The  results  indicated  that  the  derived  ANN/LES  yielded  a  predicted  pressure 
field  over  time  which  was  within  1%  of  the  experimental  data  for  all  the  a+  cases  at  all  the 
span  locations.  Further,  the  results  indicated  that  the  ANN/LES  could  accurately  extrapolate 
to  any  non-dimensional  pitch  rate  between  0.0001  and  0.2  and  to  any  span  location  from 
the  wing  root,  0%,  to  near  the  wing  tip  at  80%  span.  And,  in  all  cases,  the  linear  equation 
system  yielded  identical  results  to  those  obtained  using  the  ANN.  Thus,  it  was  possible  to 
mathematically  quantify  the  unsteady  flow  fields  obtained  experimentally.  The  techniques 
described  contribute  significantly  to  the  computational  methods  available  for  modeling 
three-dimensional  unsteady  flow  fields. 
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A  NEURAL  NETWORK  MODEL 

OF  THE  UNSTEADY  AERODYNAMICS  ON  A  PITCHING  WING 

Wiliiap'  E.  Fuller 

INTRODUCTION 

A  large  number  of  studies  have  looked  at  the  unsteady  separated  flows  associated 
with  sinusoidally  or  constant-rate  pitched  airfoils  (Robinson  and  Luttges,  1983;  Adler  and 
Luttges,  1985;  Ashwonh  et  al.,  1986;  Ashworth  and  Luttges,  1986;  Robinson  and  Luttges, 
1986;  Robinson  et  al..  1986;  Flelin  et  al.,  1986;  Robinson  and  Wissler,  1988;  Schreck  and 
Luttges,  1988;  Ashwonh  et  al.,  1989;  Schreck  and  Luttges,  1989;  Huyer  et  al.,  1990; 
Klinge  et  a’.,  1990;  Horner  et  al.,  1990;  Klinge  et  al.,  1991;  Huyer  and  Luttges,  1991; 
Schreck  et  al.,  1991;  Schreck  and  Helin,  1992)  These  studies  have  characterized  the 
unsteady  aerodynamics  using  flow  visualization  techniques,  hot-wire  annemometry  and 
surface  pressure  readings.  Funher,  these  studies  have  looked  at  both  two-  and  three- 
dimensional  unsteady  flow  fields. 

Instances  of  very  high  lift  have  been  correlated  with  the  generation  and  existence  of 
a  leading-edge  vortex  on  the  upper  surface  of  the  airfoil.  The  time  history  during 
convection  of  the  leading-edge  vortex  in  turn  determines  the  amount  of  lift  and  moments 
generated.  Further,  the  surface  pressure  and  lift  distributions  on  the  airfoil  are  not  the  same 
at  the  wing  root  and  tip.  The  interaction  between  the  leading-edge  vortex  and  the  w'ing  tip 
vortex  has  been  shown  to  be  a  highly  three-dimensional  phenomenon.  While  such  time  and 
space  dependent  changes  in  both  the  vortex  dynamics  and  lift  profile  are  desirable  from  an 
applications  standpoint  they  pose  significant  difficulties  in  both  the  prediction  and  control 
of  these  transient  factors. 

One  possibility  to  overcome  these  difficulties  might  be  to  learn  the  time  and  space 
dependencies  of  leading-edge  vortex  generation  and  convection  using  an  artificial  neural 
network  (ANN).  Neural  networks  could  then  be  utilized  as  a  model  of  the  flow  field  which 
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would  function  effectively  across  a  wide  range  of  flight  regimes.  Then,  depending  on  the 
effectiveness  of  the  neural  network  in  predicting  the  unsteady  flow  field  it  might  be 
possible  to  attempt  to  control  this  phenomenon.  However,  before  any  type  of  control 
system  can  be  attempted,  it  is  first  necessary  to  determine  the  efficacy  of  neural  networks  in 
describing  and  predicting  three  dimensional  unsteady  flow  fields. 

The  use  of  neural  networks  for  both  system  identification  and  as  control  systems  is 
emerging  as  one  possible  technique  for  handling  complex  non-linear  systems.  The  use  of 
neural  networks  to  serve  as  system  models  has  been  addressed  by  (Chu  et  al.,  1990;  Chen 
et  al.,  1990;  Ljung,  1991;  Parlos  et  al.,  1991).  These  studies  have  shown  that  non-linear 
models  of  complex  systems  can  be  developed  using  ANNs.  The  use  of  neural  networks  in 
the  design  of  control  systems  has  also  been  addressed  (Nguyen  and  Widrow,  1990; 
Narendra  and  Mukhopadhyay,  1992;  Sanori  and  Antsaklis,  1992).  Using  a  neural-network 
system  model  a  second  neural  network  is  trained  to  control  the  emulator  (model).  In  this 
fashion,  a  neural  network  can  be  trained  to  solve  a  highly  non-linear  control  problem.  More 
recently  the  use  of  neural  networks  in  the  design  of  aircraft  control  systems  has  been 
addressed  (Ha,  1991;  Troudet  et  al.,  1991;  Linse  and  Stengel,  1992;  Steck  and  Rokhsaz, 
1992).  In  these  examples,  neural  networks  were  trained  to  correlate  functions  such  as  stick 
position  with  various  aerodynamic  coefficients.  Thereby,  providing  the  opportunity  to 
implement  parts  of  an  aircraft  control  system  using  neural  networks. 

The  work  described  herein  addresses  the  difficulties  associated  with  learning  the 
complex  time  and  space  dependencies  of  leading-edge  vortex  generation  and  convection 
using  a  neural  network  architecture.  As  shown,  in  the  results,  a  highly  successful  system 
model  can  be  developed  using  ANNs.  Further,  as  shown,  a  mathematical  model  of  the 
pressure  gradient  field  over  the  surface  of  the  airfoil  can  be  derived  from  the  ANN  weight 
matrices.  This  work  makes  a  significant  contribution  to  botli  the  computational  techniques 
available  for  modeling  three-dimensional,  vortex  dominated,  unsteady  flow  fields  as  well 
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as  lo  techniques  which  should  facilitate  the  control  of  this  phenomenon. 


METHODS 

The  data  acquisition  system  is  shown  schematically  in  Fig.  1.  A  straight  wmg 
having  a  NACA  0015  cross-section  and  rectangular  planform  was  attached  to  a  circular 
splitter  plate.  Starting  at  0  degrees  this  configuration  was  dynamically  pitched  at  a  constant 
rate  about  the  wing  quarter  chord  to  an  angle  of  60  degrees  which  exceeded  the  static  stall 
angle.  During  the  pitch  history  surface  pressure  readings,  in  the  form  of  pressure 
coefficients,  were  obtained  from  15  pressure  transducers  spaced  between  0  and  90^^ 
chord.  A  total  of  54  data  records  were  obtained  which  covered  6  non-dimensional  pitch 
rates  (a+)  ranging  between  0.0001  and  0.2  and  9  span  locations  ranging  between  0%,  the 
wing  root  at  the  splitter  plate,  and  80%  span,  near  the  wing  tip.  Each  data  record  was 
comprised  of  a  total  of  200  points  which  covered  the  full  duration  of  the  pitch  cycle.  As 
shown  in  Fig.  1,  for  each  data  sample  acquired,  all  15  pressure  pons  readings  were 
simultaneously  stored  as  pressure  coefficients.  The  recorded  spatial  and  temporal  histories 
of  the  unsteady,  vortex  dominated  flows  were  then  used  to  develop  an  anificial  neural 
network  (ANN)  model  of  the  unsteady  flow  field.  The  ANN  model  was  then  used  to 
mathematically  quantify  the  three-dimensional,  vonex  dominated,  unsteady  aerodynamics 
of  a  NACA  0015  airfoil  pitched  at  constant  rates  beyond  the  static  stall  angle. 

To  model  the  unsteady  flow  field  a  paradigm  based  on  the  backpropagation  learning 
algorithm  was  developed.  This  is  schematically  shown  in  Fig.  2.  The  objective,  in  this 
case,  of  the  training  paradigm  was  to  model  the  pressure  gradient  field  using  an  ANN. 
Thus,  post-training  this  model  (ANN)  could  be  used  to  predict  the  pressure  coefficients  at 
time  (t-fA)  given  the  pressure  coefficients  at  any  time  (t).  In  general,  aerodynamic 
parameters  of  interest  which  can  be  characterized  in  this  fashion  would  include  the  time- 
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varying  surface  pressure  profile  (pressure  coefficients),  the  integrated  aerodynamic 
parameters  (lift,  drag  &  moments)  as  well  as  leading-edge  vonex  initiation,  convection  and 
shedding. 

As  shown  schematically  in  Fig.  2  a  feed-forward  architecture  with  2  hidden  layers 
was  used  which  had  the  following  configuration.  The  ANN  had  15  inputs  for  the  pressure 
coefficients  on  the  upper  surface  of  the  airfoil  (Cpi  -  Cpis).  Each  hidden  layer  was 
comprised  of  32  units  and  the  output  layer  was  comprised  of  15  units.  Bias  units  were 
included  for  each  of  the  two  hidden  layers.  A  "pattern  association"  paradigm  was  used 
where  the  input  to  the  network  at  time  (t)  was  used  to  predict  the  output  at  time  (t+A).  In 
this  case,  the  input  was  the  15  pressure  coefficients  at  time  (t)  and  the  targeted  output  was 
the  15  pressure  coefficients  at  time  (t-t-A).  For  each  pitch  history  199  consecutive  pressure 
changes  were  required  to  be  "learned"  by  the  ANN  over  the  full  pitch  cycle.  Thus,  the 
ANN  had  to  "learn"  the  pressure  gradient  field  necessary  to  generate  the  time-varying 
pressure  profiles  recorded  experimentally.  Training  was  based  on  a  supervised  gradient 
descent  algorithm,  baclqiropagation,  where  the  training  set  was  comprised  of  5  data  sets,  (5 
non-dimensional  pitch  rates  at  a  location  40%  span  from  the  splitter  plate).  The  learning  rate 
was  11=0.05  for  all  layers,  momentum  was  a=0.2  for  both  hidden  layers  and  a=0.0  for  the 
output  layer.  During  training  the  5  data  sets  were  presented  randomly  with  the  stipulation 
that  each  data  set  be  presented  an  equal  number  of  times.  The  initial  weights  were  set 
randomly  between  -0.25  and  0.25  and  training  was  performed  until  the  sum-squared  error 
was  less  than  1%  for  all  the  training  sets. 

A  novel  approach  was  taken  in  defining  the  activation  functions  utilized  in  the 
network  architecture.  A  modified  quasi-linear  function  was  used  which  had  the  following 
characteristics. 

(1)  y  =  0  and  Wij  =  0.5*Wij  forx<-a  (a  =  4.0) 

(2)  y  =  x/2a  +  0.5  for  -a  <  x  <  a  (a  =  4.0) 
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for  X  >  a 


(a  =  4,0) 


(3)  y=l  and  Wjj  =  0.5*Wij 
As  implied  by  the  equations  these  units  were  not  permitted  to  maintain  a  saturated,  y=fl  or 
y=i,  output  during  training.  This  was  implemented  by  specifying  that  all  weighted  inputs 
(Wjj)  to  a  saturated  unit  were  to  be  halved  (Wjj  =  0.5*Wjj).  Thus,  following  trainitig  all 

activation  values,  for  all  units,  were  within  the  linear  region  of  the  activation  function,  Eqn. 
2.  Since  the  output  of  each  unit  was  linear,  the  contribution  of  any  input  unit  to  any  output 
unit  remains  linearly  separable  from  all  other  inputs.  Therefore,  the  contribution  of  each 
individual  input  unit  to  each  individual  output  unit  can  be  uniquely  determined.  Thus,  post¬ 
training,  it  was  possible  to  determine,  from  the  weight  matrices  of  the  ANN,  a  single 
coefficient  which  described  the  contribution  of  each  input  unit  to  each  output  unit.  In  other 
words,  the  weight  matrices  could  be  compressed  into  a  coefficient  matrix  fA)  plus  a 
constant  vector  [K], 

(4)  [A]  [Cp(t)]  +  [k]  =  [Cp(t-hA)] 

This  in  turn  is  nothing  more  than  a  linear  equation  system  (LES)  which  acts  upon  the  input 
vector,  the  pressure  profile  at  time  (t),  imposes  the  derived  pressure  gradient  field  [A]  and 
yields  the  pressure  profile  at  time  (t-f-A).  Thus,  the  ANN  model  facilitates  the  capability  to 
mathematically  quantify  the  unsteady  aerodynamics  of  a  NACA  0015  airfoil  pitched  at 
constant  rates  beyond  the  static  stall  angle.  Since,  based  on  the  technique  employed,  the 
ANN  and  LES  must  yield  absolutely  identical  solutions,  during  all  funher  analyses  the 
ANN/LES  were  used  interchangeably  to  predict  the  time-varying  pressure  profiles. 

Post-training  to  evaluate  the  performance  of  the  ANN/LES  each  of  the  54  data 
records  was  tested.  Sum-squared  errors  were  calculated  for  each  individual  unit.  These 
values  were  then  summed  and  averaged  to  yield  a  mean  sum-squared  error  per  unit.  To 
determine  the  time-varying  error,  linear  correlations  between  the  predicted  pressure  profiles 
and  the  experimental  data  were  calculated  for  each  output  unit.  Again,  a  mean  value  was 
calculated  for  each  unit  yielding  a  mean  correlation  (r).  And,  the  performance  of  the 
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ANN/LES  was  verified  graphically  by  co-plotting  the  recorded  pressure  profiles  (raw 
data),  the  predicted  pressure  profiles  and  the  difference  between  the  tw'o  plots  (raw  data  - 
predicted).  Thus,  post-training,  it  was  possible  to  determine  not  only  how  well  the 
ANN/LES  predicted  the  training  set  data,  but  how  well  the  ANN/LES  could  extrapolate 
(generalize)  to  both  other  non-dimensional  pitch  rates  and  to  other  span  locations  not  used 
in  constructing  the  model. 

RESULTS 

Post-training  the  weight  matrices  of  the  ANN  were  used  to  derive  a  linear  equation 
system  (LES)  of  the  type  shown  in  Eqn.  4.  The  equation  system  is  shown  in  Fig.  3.  As 
shown  the  LES  acts  upon  an  input  vector,  the  pressure  profile  at  any  time  (t),  imposes  the 
derived  pressure  gradient  field  [A]  and  yields  the  pressure  profile  at  the  corresponding  time 
(t-l-A).  An  explicit  relationship  has  been  derived  where  the  influence  of  the  pressure  gradient 
field  on  each  pressure  coefficient  can  be  described  by  the  constant  coefficients  in  the  [A] 
matrix.  Thus,  post-training  either  the  ANN  model  or  the  LES  model  can  be  used  to  predict 
the  pressure  coefficients  at  time  (t+A)  given  the  pressure  coefficients  at  any  time  (t).  Note: 
since  the  ANN  and  LES  yield  absolutely  identical  solutions  in  all  cases,  during  all  further 
analyses  the  ANN/LES  were  used  interchangeably  to  predict  the  time-varying  pressure 
profiles. 

To  test  whether  or  not  the  artificial  neural  network  {ANN)/linear  equation  system 
(LES)  accurately  described  both  the  generation  and  convection  of  the  leading-edge  vonex 
as  well  as  the  highly  three-dimension  nature  of  the  flow  field,  the  equation  system  was 
tested  against  all  the  available  data  sets.  The  graphical  analysis  for  predicting  a  non- 
dimensional  pitch  rate  of  0.01  at  the  40%  span  location  are  shown  in  Figs.  4,5  and  6.  In 
each  figure  the  upper  left  hand  corner  depicts  the  pitch  history  of  the  airfoil.  The 
instantaneous  angle  of  attack,  in  degrees,  is  shown  along  the  ordinate  and  non-dimensional 
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lime  is  along  the  abscissa.  Each  figure  is  labeled  with  the  pressure  port  number  I'rom  which 


the  data  was  obtained.  Port  1  is  at  the  leading  edge  of  the  airfoil  and  pen  15  is  at  90'76 
chord,  near  the  trailing  edge  of  the  airfoil.  In  each  figure  the  ordinate  is  the  pressure 
coefficient  and  the  abscissa  is  non-dimensional  time.  In  all  figures  die  e.xperimental  (raw) 
data  is  shown  as  a  solid  line,  the  predicted  data  as  a  dotted  line  and  the  error  (raw  - 
predicted)  as  a  solid  line.  Figure  4  shows  the  results  for  the  prediction  of  the  pres.sure 
coefficients  obtained  from  pons  1-5.  Figure  5  shows  the  results  for  the  prediction  of  the 
pressure  coefficients  obtained  from  pons  6-10.  And,  figure  6  shows  the  resuits  for  the 
prediction  of  the  pressure  coefficients  obtained  from  pons  11-15.  In  ail  cases,  based  on  the 
sum-squared  error,  the  equation  system  yielded  a  predicted  pressure  field  over  time  which 
was  within  1%  of  the  experimental  data.  Table  1. 


Table  1.  Shown  is  the  sum-squared  error  for  the  predicted  pressure  coefficients  at  each 
pressure  port  as  well  as  the  conelational  coefficient  (r)  for  a  linear  correlation  between  the 
experimental  and  predicted  pressure  profiles. 


port 

ssq 

correlatio.i  (r) 

1 

1.91E-03 

0.996 

3 

3.53E-03 

0.994 

5 

2.33E-03 

0.983 

7 

1.03E-03 

0.970 

9 

1.24E-03 

0.939 

11 

6.64E-04 

0.978 

13 

7.28E-04 

0.986 

15 

6.56E-04 

0.993 

pert 

ssq 

correlation  (r) 

2 

5.19E-03 

0.993 

4 

4.03E-03 

0.988 

6 

2.63E-03 

0.962 

8 

1.31E-03 

0.940 

10 

7.53E-04 

0.965 

12 

l.OlE-03 

0.975 

14 

3.59E-04 

0.995 

The  average  correlation  (r)  was  0.977±  i.87E-02  and  the  average  sum-squared  error  for 
each  unit  was  1.82E-03±  1.44E-03.  Clearly,  the  model  can  accurately  predict  the  time 
varying  pressure  profiles. 

Moreover,  even  though  the  ANN  had  only  been  trained  on  a  subset  of  the  available 
a+  cases  at  one  span  location,  the  results  indicated  that  the  ANN/LES  was  accurate  for  all 
the  a+  cases  at  all  the  span  locations.  Funher,  the  results  indicated  that  ihe  equation  system 
could  accurately  extrapolate  to  any  non-dimensional  pitch  rate  lietwecn  0.0001  and  0.2  and 
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to  any  span  locaiion  from  ihe  wing  root,  0%,  to  near  the  wing  up  at  SO'Tc  span.  Con.sisicnt 
results  w'ere  obtained  for  ANNs  u-ained  at  different  span  locations.  In  addition,  lo  funJier 
quantify  these  results  a  4-layer  non-linear  ANN,  based  on  a  standard  sigmoidal  activation 
function,  was  trained  in  an  identical  fashion.  The  resuits  from  the  non-linear  model  were 
then  compared  with  those  of  the  linear  equation  system.  While  increased  performance  was 
noted  using  a  non-linear  network,  overall  the  pciiormance  as  measured  by  the  sum-squared 
error  increased  by  roughly  0.5%.  Thus,  the  linear  equation  system  (LES)/  ANN  was 
shown  to  provide  a  model  which  was  highly  accurate  and  showed  no  appreciable  difference 
in  performance  as  compared  to  a  non-linear  model  of  the  unsteady  flow  field. 

DISCUSSION 

A  novel  technique  was  developed  and  tested  for  computationally  modeling  three 
dimensional,  vortex  dominated,  unsteady  flow  fields  using  anificial  neural  networks 
(ANN).  Using  backpropagation  and  a  modified  quasi-Iinear  activation  function  an  ANN 
model  of  the  unsteady  flow  field  for  a  NACA  0015  airfoil  pitched  at  constant  rate  was 
derived.  The  ANN  model  was  trained  to  "learn"  the  pressure  gradient  field  underlying  the 
time-varying  pressure  profiles  recorded  experimentally.  The  ANN  model  was  then  used  to 
mathematically  quantify  the  three-dimensional,  vortex  dominated,  unsteady  aerodynamics 
of  the  phenomenon.  A  linear  equation  system  (LES)  which  yielded  identical  results  to  the 
ANN  model  was  derived  from  the  weight  matrices  of  the  ANN.  Thus,  post-training  either 
the  ANN  model  or  the  LES  model  could  be  used  to  predict  the  pressure  coefficients  at  time 
(t+A)  given  the  pressure  coefficients  at  any  time  (t). 

The  method  was  tested  on  a  total  of  54  data  records  which  covered  6  non- 
dimensional  pitch  rates  (a+)  ranging  between  0.0001  and  0.2  and  9  span  locations  ranging 
between  0%,  the  wing  root  at  the  splitter  plate,  and  80%  span,  near  the  wing  tip.  The 
results  indicated  that  the  derived  ANN  yielded  a  predicted  pressure  field  over  time  which 
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was  within  1%  of  the  experimental  data  for  all  the  a'*'  cases  at  ail  the  span  iccutions. 
Funher,  the  results  indicated  that  the  ANN  could  accurately  cxtrapoiate  lo  any  non- 
dimensional  pitch  rate  between  0.0001  and  0.2  and  to  any  span  location  from  the  wing 
root,  0%,  to  near  the  wing  tip  at  80%  span.  Funher,  in  all  cases  the  results  indicated  that 
the  LES  yielded  identical  results  to  those  obtained  using  the  ANN.  Based  on  these  results, 
it  is  reasonable  to  believe  that  not  only  can  neural  networks  be  used  to  model  unsteady  flow 
fields,  but  that  they  can  provide  new  insights  into  the  underlying  physics  of  three- 
dimensional  unsteady  aerodynamics. 

Based  on  these  results  it  is  hypothesized  that  the  underlying  physics  of  the  problem 
must  be  contained  in  both  the  ANN  and  linear  equation  system.  Thus,  the  modified  quasi- 
linear  ANNs  described  herein  provide  an  analytical  tool  with  which  to  mathematically 
identify  the  physics  of  three-dimensional,  vortex  dominated,  unsteady  flow  fields.  Funher, 
since  explicit  equation  systems,  which  are  unattainable  by  other  means,  can  be  derived 
using  the  paradigm  described,  the  development  of  system  models  and  control  systems 
targeted  at  controlling  the  behavior  of  the  leading-edge  vortex  on  a  pitching  wing  should 
now  be  possible.  In  addition,  this  approach  should  be  equally  applicable  to  other  types  of 
data  as  well  as  a  large  number  of  control  problems  where  it  is  either  very  difficult  or 
impossible  to  derive  a  set  of  linear-control  laws  for  the  system  being  modeled.  And,  this 
technique  should  make  it  possible  to  derive  a  linear  equation  system  which  approximates 
the  output  of  any  non-linear  artificial  neural  network. 
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Abstract 

The  microdynamics  of  EtAlCl,  containing  melts  are  examined  by  '^C  NMR  relaxation 
methods  as  a  function  of  melt  composition  and  temperature.  Application  of  the 
Dual  Spin  Probe  (DSP)  method  to  these  systems  reveals  interaction  between  (1)  the 
MEI*  methyl  group,  (2)  the  terminal  CHj  of  the  MEI"^  ethyl  group,  and  various 
EtAlCl^  containing  species.  Unlike  MEICl-AlClj  room  temperature  melts,  there  is 
no  indication  of  interaction  between  the  MEI'*^  ring  CH's  and  EtAlCl,. 
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INTRODUCTION 


Room  temperature  molten  salts  consisting  of  mixtures  of  AlCl,  and  l-ethyl-3- 
methylimidazolium  chloride  (MEICl),  are  of  interest  as  aprotic  solvents  for 
studying  a  wide  range  of  both  organic  and  inorganic  compounds  (1-7].  These 
chloroaluminate  molten  salts  possess  considerable  potential  as  battery 
electrolytes  and  various  types  of  electrochemical  agents  (8-10). 

The  composition  of  a  chloroaluminate  melt  has  a  considerable  effect  on  its 
physical  properties.  The  variations  in  physical  properties  of  the  melt  are  due 
to  a  combination  of  factors  including  ion- ion  interactions  (4),  and  Lewis  acid- 
base  properties.  Chloroaluminate  melts  with  AlCl,  present  in  excess  (mole 
fraction,  N,  of  AlCl,  >  0.5)  are  termed  acidic  with  AlCli  and  Al;Cl7‘  the 
predominant  anions. 

The  use  of  NMR  relaxation  methods  provides  useful  information  about  the 
dynamics  and  structure  of  various  chemical  systems  and  chloroaluminate  systems 
in  particular.  In  a  previous  workfll],  '^C  NMR  relaxation  measurements  were  used 
to  investigate  the  motion  and  interactions  of  the  MEI  cation.  The  results 
indicate  that  AlCl/  in  a  Na^^aMEI VtjAICI,  melt  forms  a  complex  by  interacting 
with  the  C-2,  C-4  and  c-5  hydrogens  on  the  MEI"^  ring.  This  investigation  was 
followed  by  studies  (12,13]  in  which  the  Dual  Spin  Probe  method  [14]  supported 
the  existence  of  MEI  (A1C1J„‘°‘'’'  complexes  in  neutral  (AlClj  =  MEICl)  and  NaCl- 
buffered  melts.  ^^Al,  “Na  and  '^C  NMR  relaxation  results  confirmed  the  presence 
of  the  chloroaluminate-MEI*  complexes  and  yielded  ^A1  and  °Na  liquid 
state  quadrupole  coupling  constants( 12, 13 J . 

Application  of  the  Dual  Spin  Probe(DSP)  relaxation  method  typically  requires 
knowledge  of  ’’c  dipolar  relaxation  rates  which  are  defined  by  (1),  the  basic 
equation  in  which  the  ’’c  nucleus  is  relaxed  by  'H[15]: 

where  R,"  (=  l/T")  is  the  dipolar  relaxation  rate,  N„  is  the  number  of 
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hydrogens  attached  directly  to  the  carbon  atom,  Yc  and  Ym  are  gyromagnetic  ratios 


and  Ten  =  1.09  A.  is  the  effective  correlation  time  and  varies  exponentially 
with  temperature.  Equation  (1)  is  operative  while  under  the  "extreme  narrowing 
condition"  (uT^(f«l)  which  is  usually  applicable  for  small  molecules  including 
the  chloroaluminate  melts(ll). 

is  ootained  by  measuring  T,,  the  Nuclear  Overhauser  Enhancement  factor, 
h  (»1m«  =  Yii/2Yc)  and  using  eqn  (2)  [16]; 

R,‘“  =  nR,/l-988  (2) 

The  other  part  of  the  DSP  method  requires  knowledge  of  quadrupoiar 
relaxation  rates  for  nuclei  such  as  ’’Al  and  ^Na.  If  there  is  a  distortion  from 
tetrahedral  or  cubic  symmetry,  nuclei  such  as  ^A1  and  ^Na  will  be  under  the 
influence  of  an  electric  field  gradient  hich  produces  the  quadrupole 
interaction.  The  quadrupoiar  relaxation  rate  in  the  "extreme  narrowing  region" 
is  given  by(3)  (15,17j; 


R,  =  (37r2(2I+3)/10I^(2I-l))[l+(2V3))(e=Qq/h)^T,  (3) 

where  I  *  3/2  for  ^Na  and  5/2  for  ^Al,  eQ  is  the  nuclear  quadrupole  moment,  eg 
is  the  maximum  component  of  the  electric  field  gradient  tensor,  and  z  is  the 
asymmetry  parameter  of  the  electric  field  gradient  tensor (z  =  0  for  AlClj) . 

The  quadrupole  coupling  constant,  QCC,  is  given  by; 

QCC  =  [e^q/h]  (4) 

The  DSP  method  has  been  applied  to  chloroaluminate  melts [12, 13]  and  has 
provided  evidence  that  the  ring  hydrogens  of  MEI'^  interact  with  the 
tetrachloroaluminate  anion.  The  existence  of  these  complexes  has  been  supported 
fay  linear  plots  of  ‘^C  dipolar  relaxation  ratesiRj"*^)  vs.  quadrupoiar  ‘^Al 
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relaxation  rates(R|)  that  pass  through  the  origin  as  predicted  by  equation  (5): 


R,“'(‘3C)/N„(>>YcY..)'rcH*  =  R,(-Al)/ax'  (5) 

where  a  =  [3;rVlO][(2I  +  -  1>]I1  -*•  (z‘/2)],  and  QCC  =  x- 

During  this  summer  research  program,  the  DSP  method  was  applied  to  melts 
containing  MEICl,  AlClj  and  EtAlCl,.  The  inclusion  of  EtAlCl,  provided  a 
"baseline"  as  there  is  a  covalent  bond  between  the  ethyl  group  and  aluminum  in 
BtAlCl;.  The  existence  of  covalent  banding(or  ccmplexation)  between  quadrupolar 
and  dipolar  nuclei  in  a  molecule  results  in  a  linear  plot  of  eqn.  (5)  that  passes 
through  the  origin.  In  the  MEICl-EtAlCl,  melts  reported  herein,  we  observed  a 
linear  plot  of  eqn  (5)  that  passed  through  the  origin  when  applied  to  the 
terminal  CH3  carbon  in  EtAlCl,  and  one  of  the  peaks  in  the  *^A1  NMR  of  the  melts. 

EXPERIMENTAL 

Materials 


The  l-ethyl-3-methylimidazolium  chloride  <MEIC1)  and  chloro-aluminate  molten 
salts  were  prepared  as  described  previously  [1].  Ethylaluminum  dichloride 
(EtAlClj)  was  obtained  from  Aldrich.  All  materials  were  stored  under  anhydrous 
helium  gas  atmosphere  in  a  dry  box.  All  molten  salt  preparations  and 
manipulations  were  performed  in  the  dry  box.  Samples  were  loaded  into  5  mm 
sample  tubes,  capped  in  the  dry  box,  removed,  and  sealed  immediately  with  a 
torch. 

NMR  Measurements 

'^C  NMR  spin-lattice  relaxation  times  were  recorded  this  summer  on  a  Varian 
Gemini-300  spectrometer  at  75.43  MHz  and  ^Al  NMR  spin-lattice  relaxation  times 
were  recorded  previously  on  a  Varian  XL-300  spectrometer  operating  at  78.15  MHz. 
Temperature  measurements  were  calibrated  against  methanol  or  ethylene  glycol  and 
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are  accurate  to  within  O-S^C.  Pulse  widths(90°)  were  typically  8.6  {75.43  KHz) 
and  7.6(78.15  MHz)  ps.  Longitudinal  relaxation  times  were  measured  by  the  the 
inversion-recovery  method  ( 180°-t-90°-T)  with  T>10T,.  At  least  12  delay  times(T) 
were  used  and  the  results  fitted  to  a  three  parameter  exponential.  ‘’c  NOE 
measurements  were  made  using  the  gated  decoupler  method(18).  It  is  likely  that 
the  error  in  the  NOE  measurements  is  in  the  5-10%  range(18]. 

RESULTS  AND  DISCUSSION 

The  ability  of  both  AlClj  and  EtAlCl,  to  form  dimers  [  19 , 20  ]  led  us  to 
examine  the  ’’Al  spectra  of:  (1)  neat  EtAlCl;,  (2)  MEICl-EtAlCl,  and  (3)  ternary 
melts  (N  =  AlClj/MEICi/EtAlCl,)  [21  ] .  The  neat  EtAlCl,  ’’Al  NMR  spectrum  contains 
two  peaks  (21).  Peak  1  is  a  broad  downfield  peak  that  dominates  the  spectrum. 
The  second  peak  (upfield)  overlaps  peak  1  and  is  only  a  fraction  of  peak  1  in 
total  peak  area.  Peak  2  collapses  into  peak  1  as  the  temperature  is  lowered  from 
60  to  25 ®C.  These  two  aluminum  sites  are  consistent  with  the  extent  of  monomer- 
dimer  formation  in  liquid  EtAlCl2(21). 

The  MEICl-EtAlClj  (N  =  0.5/0. 5)  melt  ^A1  NMR  spectrum  also  has  two  peaks. 
In  this  case,  peak  1  (downfield)  is  very  broad  while  peak  2  is  very  sharp,  and  has 
a  low  peak  area.  Peak  2  increases  slightly  in  area  and  peak  1  broadens  as  the 
temperature  is  lowered  from  70  to  O^C.  We  have  previously ( 21 )  made  the  tentative 
assign-ments  of  EtAlClj"  for  peak  1  (downfield)  and  EtjAl^Clj'  for  peak  2. 

In  this  study,  we  first  apply  the  DSP  method  to  the  CH,  carbon  in  EtAlClj  and 
^Al  NMR  peaks  1  and  2  from  several  melt  combinations  and  neat  EtAlClj.  Fig.  1 
contains  the  results  for  ^A1  peak  l(downfield)  and  Fig.  2  contains  the  results 
for  ^Al  peak  2.  The  fact  that  both  plots  are  linear  and  pass  through  the  origin, 
indicate  that;  (1)  the  DSP  method  is  appropriate  for  these  systems  and  (2)  the 
species  associated  with  each  peak  contains  EtAlCl,.  Furthermore,  the  slopes  of 
these  lines  can  be  used  to  calculate  the  relative  quadrupole  coupling  constants 
for  the  EtAlClj-containing  species  in  solution. 
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Fig.  1.  Dipolar  Rl's  vs  '^Al  Rl's(25  to  TO'C)  for  A1  peak  1  (127-131  ppm  from 
A1(H,0)/")  . 
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Fig.  2.  "C  Dipolar  Rl's  vs  ^A1  Rl'8<25  to  70*C)  for  A1  peak  2  (102.5-103.0  ppn 
from  A1(H,0)4’^)  . 


The  QCC  values  obtained  from  Fig.  1(A1  peak  1)  are  171,  119,  106  and  93  MHz 
for  the  (.5/.  5),  ( . 35/ . 40/ . 25  )  ,  ( . 2 5/ . 40/ . 3 5  )  melts  and  neat  EtAlCl., 
respectively.  The  QCC  values  obtained  from  Fig.  2(A1  peak  2)  are  6.9,  20,  11  and 
93  MHz  for  the  (.5/. 5),  {  ,  35/ . 40/ . 25 ) ,  (  .  2 5/ . 40/ . 35 )  melts  and  neat 
EtAlCl.(repeated)  . 

Results  of  the  Dual  Spin  Probe  method  (eqn.  [5])  applied  to  the  (.5/. 5), 
( .  35/ . 40/ . 25 )  and  (  .  2 5/ . 40/ . 35 )  melts  indicate  interactions  between  the  Al- 
containing  species  in  peak  2  { 102 . 5-103 . 0  ppm  relative  to  Al(H,0)j’')  and  both  the 
NCH,  and  ethyl  terminal  CH,  groups  of  Fig.  3  contains  the  plots  for  the 
NCH,  group  in  each  melt  and  Fig.  4  contains  data  for  the  terminal  CH,  on  the  MEI 
ethyl  group. 

The  QCC's  obtained  from  the  slopes  in  Fig.  3(HEI  NCH,)  are  1.7,  2.3  and  4.4 
MHz  for  the  (.5/. 5),  ( . 35/ . 40/ . 25 )  and  ( . 25/ . 40/ . 35)  melts.  The  QCC's  from  Fig. 
4{terminal  CH,  on  the  MEI  ethyl  group)  are  1.6,  6.9  and  1.3  MHz  for  the  (.5/. 5), 
( .35/. 40/. 25)  and  ( . 2 5/ . 40/ . 35 )  melts. 

Finally,  there  is  no  correlation  between  the  ring  hydrogen  dipolar  Rl's  and 
any  of  the  ''Al  peak  Rl's.  This  result  is  directly  opposite  to  that  found  in 
MEICl-AlCl,  systems  [11,12]. 

CONCLUSIONS 

Application  of  the  DSP  method  to  these  mixed  melt  systems  indicates  a  lack 
of  complexation  between  the  ring  hydrogens  of  MEI*  and  any  of  these  aluminum 
containing  species.  These  and  previous  results [21]  suggest  that  the  formation 
of  various  charged  dimers  containing  EtAlClj  takes  precedence  over  the  formation 
of  complexes  between  EtAlCl,"  and  the  MEI*  ring  hydrogens.  However,  there  is 
evidence  of  interaction  between  the  various  Al-containing  species  and  the  CH, 
groups(NCH,  and  terminal  CH,  in  the  ethyl  group)  of  MEI*  in  these  melts. 
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Abstract 

The  aim  of  this  study  was  to  develop  both  a  laboratory  model  of  closed  head-injury  and  an  analytical 
model  of  venous  blood  flow  from  the  brain  to  test  the  hypothesis  that  variations  in  venous  pressure  associated 
with  the  respiratory  cycle  can  have  a  dominant  influence  on  venous  flaw  from  the  brain  during  elevated 
intracranial  pressure.  A  young  adult  pig  with  an  implanted  intraaaniai  balloon  designed  to  manipuiate 
intracranial  volume  was  used  as  a  laboratory  mod^.  Art  analog  ^ectricai  circuit  model  was  used  to  provide 
a  theoretical  analytical  description  of  cerebral  venous  blood  flow  during  elevated  intracranisfl  pressure. 
Both  experimerrtal  and  theoretical  results  indicate  that  during  intact  autoregulation  of  cerebral  blood  flow, 
respiratory  induced  venous  pressure  changes  systematically  influence  intracranial  blood  volume. 
Specifically,  intracranial  blood  volume  increases  during  inhalation  and  decreases  during  expiration. 
Furthermore,  the  difference  in  change  of  intracranial  volume  between  the  two  phases  of  ventilation, 
inhalation  and  expiration,  increases  with  increasing  mean  intracranial  pressure.  However,  during  loss  of 
regUation  of  cerebral  blood  flow,  venous  blood  flow  and  the  resulting  chan^  of  intracranial  blood  volume 
are  not  systematically  influerKed  by  respiratory  induced  verxxjs  pressure  changes. 


